1/2 


AD-A145  382  A  FINITE  ELEMENT  PROJECTION  METHOD  FOR  THE  SOLUTION  OF 
PARTICLE  TRANSPORT.  .  (U)  AIR  FORCE  INST  OF  TECH 
MRIGHT-PATTERSON  AFB  OH  E  E  HILLS  JUL  84 
UNCLASSIFIED  AFIT/CI/NR-84-56D  F/G  12/1  NL 


nor  rf)pY  AD- A 145  382 


SECURITY  CLASSIFICATION  or  THIS  PACE  (Whtn  Dtf'Enertd)' 

REPORT  DOCUMENTATION  PAGE  bworecompletSo tor.. 

i.  “report  NUMBER  -  |2.  GOVT  ACCESSION  NO.  3.  RECiPlEMX'S  CATALOG  NUMBER 


i.  REPORT  NUMBER  2.  GOVT  ACCESSION  NO.  i.  RECJrifAT  5  CATALOG  NUMBER 

AFIT/CI/NR  84-56D  A» -  A  S 

«.  TITLE  (m>d  Submit)  *•  TYPE  OF  REPORT  ft  PERIOD  COVEREO 

A  Finite  Element  Projection  Method  For  The  THESIS/DISSERTATION 

Solution  Of  Particle  Transport  Problems  With  mi/l _ 

Anisotropic  Scattering  performing  org.  report  number 


1 7.  AUTHORS 


Eze  Ewart  Wills 


contract  or  grant  number^ 


9.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

AFIT  STUDENT  AT:  The  University  of  New  Mexico 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  ft  PORK  UNIT  NUMBERS 


July  1984 


IS.  NUMBER  OF  PAGES 


II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  «*•  REPORT  DATE 

AFIT/NR  July  1984 _ 

WPAFB  OH  45433  '»•  OF  PACES 

14.  MONITORING  AGENCY  NAME  4  AOORES S(ll  dllltranl  Irom  Controlling  Ollict)  IS.  SECURITY  CLASS,  (ol  thit  report; 

UNCLASS 

15*.  DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


[16.  DISTRIBUTION  STATEMENT  (at  tbit  Report; 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


ELECTE| 
SEP  1  3  1984 


'17. '  OlS'OtlBUTION  STATEMENT  (ol  the  ebetract  entered  In  Block  20,  II  dllimrent  from  Report) 


lie.  SUPr’l  LMENTARY  NOTES 


i/VWVH)  FOR  PUBLIC  RELEASE:  IAW  AFR  190-1 


LY$N  E.  WOLAVER 
Dean  for  Research  and 
S' 4 Professional  Developmeni 
_ AFIT,  Wright-Patterson  AFB  OH 


[  19.  KEY  WORDS  (Continue  on  reverie  side  it  neceesary  and  identify  by  block  number) 


j  20.  ABS't  RAC T  ( Continue  on  reverse  aide  II  neceaaary  and  Identity  by  block  number) 


ATTACHED 


rpqu  m  ...a 

i  j  1473  EDITION  OF  1  NOV  65  IS  OBSOLETE 

84  09  13  021 


UNCLASS 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  Dttt  Erttrri \ 


A  FINITE  ELEMENT  PROJECTION  METHOD 
FOR  THE  SOLUTION  OF  PARTICLE  TRANSPORT 
PROBLEMS  WITH  ANISOTROPIC  SCATTERING 


EZE  EWART  WILLS 

B.S.,  University  of  Maryland,  1978 
M.S.,  Air  Force  Institute  of  Technology,  1981 


DISSERTATION 

Submitted  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 

* 

Doctor  of  Philosophy  in  Engineering 


The  University  of  New  Mexico 
Albuquerque,  New  Mexico 

July,  1984 


84  09  13  021 


AF IT/C I/NR  84-56D 


AFIT  RESEARCH  ASSESSMENT 

The  purpose  of  this  questionnaire  Is  to  ascertain  the  value  and/or  contribution  of  research 
accomplished  by  students  or  faculty  of  the  Air  Force  Institute  of  Technology  (AU).  It  would  be 
greatly  appreciated  If  you  would  complete  the  following  questionnaire  and  return  It  to: 

AF IT/NR 

Wrlght-Patterson  AFB  OH  45433 

RESEARCH  TITLE:  _ A_Finite  Element  Projection  Method  For  The  Solution  Of  Particle  Transport 

Problems  With  Anisotropic  Scattering _ 

AUTHOR:  Fzp  Ewart.  Mills _ 

RESEARCH  ASSESSMENT  QUESTIONS: 

1.  Did  this  research  contribute  to  a  current  Air  Force  project? 

(  )  a.  YES  {  )  b.  NO 

2.  Do  you  believe  this  research  topic  is  significant  enough  that  it  would  have  been  researched 
(or  contracted)  by  your  organization  or  another  agency  if  AFIT  had  not?' 

(  )  a.  YES  (  )  b.  NO 

3.  The  benefits  of  AFIT  research  can  often  be  expressed  by  the  equivalent  value  that  your 
agency  achieved/received  by  virtue  of  AFIT  performing  the  research.  Can  you  estimate  what  this 
research  would  have  cost  if  it  had  been  accomplished  under  contract  or  if  it  had  been  done  in-house 
in  terms  of  manpower  and/or  dollars? 

(  )  a.  MAN-YEARS  _  (  )  b.  $ _ 

4.  Often  it  is  not  possible  to  attach  equivalent  dollar  values  to  research,  although  the 
results  of  the  research  may,  in  fact,  be  important.  Whether  or  not  you  were  able  to  establish  an 
equivalent  value  for  this  research  (3.  above),  what  is  your  estimate  of  its  significance? 

(  )  a.  HIGHLY  (  )  b.  SIGNIFICANT  (  )  c.  SLIGHTLY  (  )  d.  OF  NO 

SIGNIFICANT  SIGNIFICANT  SIGNIFICANCE 

5.  AFIT  welcomes  any  further  corrments  you  may  have  on  the  above  questions,  or  any  additional 
details  concerning  the  current  application,  future  potential,  or  other  value  of  this  research. 

Please  use  the  bottom  part  of  this  questionnaire  for  your  statement (s). 


ACKNOWLEDGEMENT 


I  wish  to  express  thanks  to  members  of  my  dissertation 
committee  for  their  support  and  assistance  throughout  this 
research  project.  My  principal  advisors  Dr.  N.  Roderick 
and  Dr.  P.  McDaniel  provided  invaluable  guidance  and  many 
hours  of  fruitful  discussion.  Without  their  support  this 
research  would  not  have  been  possible. 

Also,  I  would  like  to  thank  the  faculty  and  staff  of 
the  Nuclear  and  Chemical  Engineering  Department  and  also 
the  Mathematics  Department  of  the  University  of  New  Mexico 
for  their  help  and  suggestions.  Special  thanks  must  be 
given  to  Dr.  R.  Allen,  Dr.  S.  Pruess  and  Dr.  B.  Russell  of 
the  Mathematics  Department  of  the  University  of  New  Mexico 
for  many  hours  of  personal  and  classroom  advice  and  dis¬ 
cussion. 

I  am  also  grateful  to  the  United  States  Air  Force  for 
providing  me  the  opportunity  to  complete  this  work. 
Specifically,  I  would  like  to  express  my  appreciation  and 
thanks  to  the  Air  Force  Weapons  Laboratory  for  their  sup¬ 
port  in  providing  access  and  use  of  their  facilities,  and 
also  for  the  encouragement,  assistance  and  advice  which 


they  provided.  I  must  also  acknowledge  and  thank  the 
faculty  and  staff  of  the  Air  Force  Institute  of  Technology 
for  recommending  this  research  area,  providing  initial 
support  and  also  convincing  me  of  the  need,  importance 
and  possibility  for  significant  research  in  this  area. 

Many  others  too  numerous  to  mention  herein  have 
contributed  to  and  supported  me  in  this  work.  To  all  of 
them  I  would  like  to  extend  a  heartfelt  and  sincere  thank 
you  and  apologize  for  not  being  able  to  individually 
acknowledge  their  support  and  assistance. 

Finally,  I  must  thank  my  family  for  their  patience 
and  moral  support.  The  sacrifice  and  understanding  of  my 
wife  and  children  were  essential  to  the  successful 
completion  of  this  work.  I  must  also  thank  my  wife  for 
typing  this  manuscript. 


Accession  For 
IJTIS  GRA&I  v 

DXIC  TAB  □ 

Unannounced  □ 

Justification - 


Distribution/ 
Availability  Codes 


A  FINITE  ELEMENT  PROJECTION  METHOD 
FOR  THE  SOLUTION  OF  PARTICLE  TRANSPORT 
PROBLEMS  WITH  ANISOTROPIC  SCATTERING 


Eze  Ewart  Wills 


ABSTRACT  OF  DISSERTATION 

Submitted  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 


Doctor  of  Philosophy  in  Engineering 

The  University  of  New  Mexico 
Albuquerque,  New  Mexico 

July,  1984 


A  FINITE  ELEMENT  PROJECTION  METHOD 
FOR  THE  SOLUTION  OF  PARTICLE  TRANSPORT 
PROBLEMS  WITH  ANISOTROPIC  SCATTERING 


Eze  Ewart  Wills 

B.S.,  Engineering,  University  of  Maryland,  1978 
M.S.,  Engineering,  Air  Force  Institute  of  Technology,  1981 
Ph.D. ,  Engineering,  University  of  New  Mexico,  1984 

"  ~^>  A  solution  method  for  solving  particle  transport 
problems  has  been  developed.  This  solution  approach 
embodies  a  finite  element  projection  technique  and  a 
related  equivalent  variational  Raleigh-Ritz  formalism. 
Particle  flux  in  the  transport  equation  is  expressed  as 
a  linear  and  separable  Siam  of  odd  and  even  components 
in  the  direction  variables.  Then  a  classical  variational 
principle  is  obtained  and  shown  to  be  equivalent  to  a 
Bubnov-Galerkin  projected  solution.  A  dual  finite  element 
basis  of  polynomial  splines  in  space  and  spherical  harmonics 
in  angle  is  used  in  the  Bubnov-Galerkin  equations . 

The  general  theoretical  and  numerical  problem  for¬ 
malism  is  carried  out  in  a  3 -dimensional  geometry  with 
anisotropic  scattering  and  with  a  piecewise  constant  energy 
dependence.  This  is  a  seven-dimensional  problem  with  time 
dependence,  three  spatial  and  two  angular  or  directional 
variables  and  with  a  multigroup  treatment  of  the  energy  — 


\ 

\ 
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''^dependence.  The  boundary  conditions  for  most  physical 
problems  of  interest  are  dealt  with  explicitly  and 
rigorously  by  a  classical  minimization  (variational)  prin¬ 
ciple.^ 

'  The  solution  method  is  developed  as  a  complementary 

alternative  to  the  standard  solution  techniques  of  Discrete 

Ordinates,  Monte  Carlo  and  the  P„  method.  The  Galerkin 

n 

projected  operator  and  transport  matrix  are  positive 
definite,  symmetric  and  self-adjoint.  This  insures 
existence,  uniqueness,  and  convergence  of  the  solution. 

This  formalism  does  not  have  the  inherent  properties  which 
\  have  produced  the  ray  effect  problem  in  discrete  ordinates. 
)  The  computational  validation  of  the  method  was 
obtained  by  a  computer  solution  to  the  air-over-ground 
problem.  This  problem  is  of  significant  interest  in  the 
areas  of  nuclear  weapons  effects  and  radiation  physics. 

It  is  modeled  in  cylindrical  (r,z)  geometry  with  an 
exponentially  varying  atmosphere,  anisotropic  scattering, 
anisotropic  first  scatter  sources,  and  with  the  air-ground 
interface  included. 

Results  for  the  air-over-ground  problem  are  presented. 
These  results  show  that  this  solution  approach  mitigates 
ray  effects.  They  also  show  the  potential  of  this 
technique  to  serve  as  a  viable  alternative  to  Discrete 
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Ordinates  and  Monte  Carlo.  However,  further  work  in 
extending  the  computer  implementation  of  the  method  to 
time  and  energy  dependent  problems,  and  to  solving  and 
validating  this  technique  on  a  larger  class  of  particle 
transport  problems  is  required. 
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CHAPTER  I 


INTRODUCTION 


Particle  transport  occurs  in  many  physical  and 
engineering  systems.  Therefore,  the  solution  of  particle 
transport  problems  is  of  widespread  engineering  and 
scientific  interest  in  many  diverse  areas  of  research, 
development  and  system  applications.  In  the  areas  of 
plasma  physics,  nuclear  energy  (fission  and  fusion); 
radiation  effects;  space  research  and  exploration,  and 
many  others,  the  theory,  physics  and  modeling  of  particle 
transport  processes  play  an  important  and  vital  role. 

There  are  numerous  solution  techniques  and  approaches 
for  solving  particle  transport  problems  (Sanchez  and 
McCormick  1981).  However,  the  most  widely  used  solution 
methods  are  based  on  a  particle  tracking  and  statistical 
treatment  of  individual  particle  interactions  and  trans¬ 
port  or  on  the  principles  of  particle  conservation  and 
distribution  functions.  In  the  latter  case  particle 
conservation  is  usually  expressed  in  terms  of  an  integro- 
differential  equation  (the  transport  equation)  and  particle 


number  densities  and  distribution  functions  are  obtained 
from  a  solution  to  this  equation. 

In  general  numerical  solutions  to  particle  transport 
problems  are  usually  obtained  by  the  Pn  method,  discrete  or¬ 
dinates  (S„) ,  or  Monte  Carlo.  The  P  and  techniques 
n  n  n 

provide  solutions  to  the  transport  equation  whereas  Monte 
Carlo  is  based  on  particle  tracking  and  statistics.  These 
methods  provide  accurate  solutions  for  a  large  number  of 
problems.  However,  for  some  problems  they  have  severe 
limitations . 

The  Monte  Carlo  method,  usually  requires  the  use  of 
many  hours  of  expensive  computer  time,  and  for  deep 
penetration  problems  it  is  subject  to  errors  due  to  statis¬ 
tical  inadequacies  (Janni,  1979).  Discrete  ordinates  may  use 
fewer  computer  resources  than  Monte  Carlo,  however  it  has 
a  computational  difficulty  called  "ray  effects"  (Lathrop, 
1968).  Ray  effects  are  a  result  of  the  angular  discreti¬ 
zation  of  the  transport  equation.  In  the  Sn  method  this 
discretization  constrains  particles  to  move  in  specific 
directions.  This,  in  turn,  produces  distortions  in  the 
particle  flux.  For  problems  with  strong  absorbers, 
localized  sources,  and  high  energy  streaming  particles, 
these  distortions  produce  solutions  which  are  inaccurate. 
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The  P^  method  or  diffusion  theory  is  the  most  often 

used  P  method.  However,  unlike  Monte  Carlo  and  the  S 
n  n 

method,  it  does  not  allow  a  detailed  treatment  of  the 
angular  dependence  of  transport  processes.  This  limita¬ 
tion  makes  diffusion  theory  inadequate  for  problems  which 
are  non-diffusive .  However,  higher  order  PR  approxima¬ 
tions  are  computationally  inefficient  and  expensive, 
especially  in  two  or  three  space  dimensions  (Kaplan,  et  al., 
1966). 

The  Problem  and  Scope 

The  problem  is  to  rigorously  develop  a  method  for 
solving  general  particle  transport  problems.  This  should 
be  a  method  that  removes  some  of  the  limitations  and 
inherent  difficulties  of  Monte  Carlo,  discrete  ordinates 
and  the  PR  method.  Furthermore,  this  solution  approach 
should  be  able  to  serve  as  a  complementary  alternative  to 
the  presently  used  techniques  for  solving  particle  trans¬ 
port  problems . 

In  order  to  accomplish  this  a  detailed  formalism  for 
solving  the  transport  equation  is  required.  Specif icially , 
a  mathematical  and  numerical  model  to  solve  the  general 
transport  equation  with  anisotropic  sources  and  scattering 
and  general  boundary  conditions  will  be  developed.  This 
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model  is  an  extension  of  current  approximate  numerical 
techniques  for  solving  particle  transport  and  other 
engineering  problems.  Also,  it  has  stable  mathematical 
and  numerical  properties. 


The  Solution  Approach 

The  solution  approach  is  to  first  develop  a  mathema¬ 
tical  formalism  (model)  by  drawing  upon  the  present  body 
of  knowledge  and  experience  in  solving  engineering  and 
physical  problems  in  general  and  particle  transport  problems 
in  particular.  We  begin  by  first  examining  current  solution 
methods  for  solving  the  transport  equation  and  then  proceed 
by  using  insight  and  hopefully  good  judgment  to  further 
develop,  examine  and  rigorously  explain  our  mathematical 
and  numerical  models . 

Our  numerical  model  is  based  on  the  well-known 
techniques  of  finite  elements,  projection  operators  and 
a  classical  Raleigh-Ritz  solution.  The  mathematical 
model  is  derived  from  a  variational  principle  and  an 
equivalent  Bubnov-Galerkin  solution  of  the  second  order 
forms  of  the  transport  equation.  These  second  order 
equations  were  first  introduced  by  Kaplan  and  Davis  in 


Our  intent  is  to  propose  a  solution  strategy  that  is 
similar  to  diffusion  theory  but  has  all  the  detailed  angu¬ 
lar  or  vector  information  of  anisotropic  scattering  and 
sources.  We  also  develop  a  rigorous  treatment  of  general 
boundary  conditions.  This  approach  differs  from  discrete 
ordinates  and  the  Pn  method  in  the  general  formalism, 
solution  strategy  and  properties  and  also  in  the  treatment 
of  boundary  conditions .  To  accomplish  this  we  take  the 
view  that  a  successful  numerical  solution  strategy 
should  be  developed  in  accordance  with  the  properties  of 
a  rigorous  mathematical  formulation.  In  fact  the  numeri¬ 
cal  and  mathematical  models  to  the  problem  should  go  hand 
in  hand.  Ideally  we  would  like  to  have  a  mathematical 
model  that  complements  and  is  compatible  with  the  numerical 
approach  and  vice  versa.  In  such  a  situation  the  numerical 
solution  strategy  when  implemented  will  be  best  suited  to 
highlight  and  use  the  special  mathematical  properties  of 
the  problem  equations  and  thus  provide  a  reasonable  or 
"best"  approximate  solution.  We  therefore  approached  the 
problem  by  first  seeking  a  mathematical  formulation  that 
had  desirable  mathematical,  computational  and  numerical 
properties.  These  properties  which  are  present  in  the 
or  diffusion  equation  are  desirable  in  that  they  would 
facilitate  a  numerical  solution  to  the  problem.  In  fact 


they  should  allow  a  solution  by  similar  and  well  established 
techniques  which  are  being  used  to  solve  diffusion  problems. 

Secondly,  we  chose  our  numerical  model  because  of  our 
mathematical  formalism  and  equations  and  the  need  to  solve 
general  problems  of  physical  and  engineering  interest.  In 
this  numerical  development  we  attempted  to  recognize  and 
use  the  properties  of  our  mathematical  model.  Finally 
we  applied  this  solution  technique  to  a  real  problem,  to 
examine,  validate  and  compare  the  results. 

In  summary,  our  mathematical  development  is  based  upon 
the  parity  and  odd-even  properties  of  spherical  harmonics 
coupled  with  a  Bubnov-Galerkin  and  classical  variational 
approach.  However,  our  numerical  development  evolves  pri¬ 
marily  upon  the  use  of  finite  elements  and  projection  opera¬ 
tors.  This  use  of  a  finite  element  projection  technique  is 
a  logical  extension  of  our  mathematical  formalism.  Further¬ 
more  it  provides  the  required  flexibility  needed  in  solving 
general  and  difficult  problems. 

Projection  Methods  and  Finite  Elements 

Over  the  last  two  decades  there  has  been  a  successful 
evolution  of  projection  methods  and  finite  elements  in 
solving  engineering  problems.  The  use  of  finite  elements 
and  in  particular  projection  methods  began  in  the  early 


19th  century.  The  finite  element  method  is  really  a  projec¬ 
tion  method  although  it  is  not  usually  recognized  as  such. 
Furthermore,  many  approximate  methods  for  solving  integral 
and  differential  equations  are  also  projection  methods.  A 
detailed  discussion  of  projection  methods  can  be  found  in 
Atkinson  (1976).  Projection  methods  include  such  widely 
used  solution  techniques  as  Galerkin,  collocation,  least 
squares  and  the  method  of  moments  (Goldberg,  1978).  A  pre¬ 
cise  mathematical  description  of  these  methods  is  presented 
in  Chapter  VI. 

Early  19th  century  mathematicians  began  using  projec¬ 
tion  methods  to  solve  ordinary  and  partial  differential 
equations,  which  could  not  be  solved  in  closed  form  (Kantor¬ 
ovich  and  Krylov,  1958).  They  developed  techniques  for 
expanding  the  solution  as  an  infinite  series  of  known 
functions  with  undetermined  coefficients.  Constraints  or 
conditions  were  then  applied  to  the  approximate  equations 
in  order  to  "determine"  the  undetermined  coefficients.  The 
use  of  different  constraints  in  determining  the  coefficients 
produced  various  projection  methods  whereas  the  use  of 
special  functions  to  expand  the  solution  resulted  in  the 
finite  element  method. 

Projection  methods  can  be  used  to  solve  linear, 
non-linear  and  non-self -adjoint  problems.  These  methods 


are  sometimes  called  the  method  of  weighted  residuals  or 
error-distribution  principles  by  engineers  and  applied 
mathematicians  (Finlayson  and  Scriven,  1966).  There  are 
many  types  and  classifications  of  projection  methods 
depending  not  only  upon  the  constraints  by  which  the 
coefficients  are  determined  but  also  on  how  the  boundary 
conditions  are  treated.  A  somewhat  unified  treatment  of 
some  of  these  methods  was  developed  by  Crandall  in  1956 
and  Finlayson  and  Scriven  in  1966.  However,  many  of 
these  approximate  solution  techniques  were  introduced  in 
the  1920’s  and  1930's.  Galerkin  introduced  his  method  as 
early  as  1915  whereas  the  concept  of  finite  elements  was 
proposed  in  1924  by  C.  B.  Biezeno.  He  recommended  the  use 
of  special  piecewise  constant  functions  in  solving  stress 
problems.  He  called  this  projection  technique  the  subdomain 
method. 

However,  it  was  in  1943  that  the  finite  element 
method  was  formally  introduced  by  Courant.  He  proposed 
the  use  of  piecewise  Lagrange  polynomials  on  a  triangular 
mesh  to  solve  two-dimensional  vibration  problems.  Subse¬ 
quently  and  with  the  advent  of  computers,  the  finite 
element  method  has  been  used  to  solve  a  wide  range  of 
problems.  One  of  the  attractive  features  of  finite 
elements  is  the  inherent  capability  to  handle  difficult 


problems  with  complicated  geometries.  General  particle 
transport  problems  of  engineering  and  physical  interest 
are  usually  in  this  category. 

Review  of  the  Literature 

The  first  use  of  a  finite  element  projection  technique 
to  solve  particle  transport  problems  occurred  in  the  early 
1970's.  In  1971  a  number  of  researchers  including  Kaper, 
et  al.,  Semenza,  et  al.f  Kang  and  Hansen,  and  Ohnishi  used 
finite  elements  to  solve  the  diffusion  and  steady  state 
transport  equation.  In  1972  a  more  detailed  examination 
of  the  use  of  finite  elements  to  solve  neutron  diffusion 
problems  was  provided  by  Kaper,  Leaf  and  Lindeman.  They 
concluded  that  a  finite  element  solution  of  the  diffusion 
equation  had  advantages  in  terms  of  accuracy  and  conver¬ 
gence  in  comparison  to  the  usual  finite  difference  solution. 

Other  researchers  including  Miller,  et  al.,  Ukai,  Hill, 
Yaun  and  Lewis,  Seed,  et  al.,  Finkelstein  and  Krumbein  and 
many  others,  have  directly  applied  finite  elements  to  solve 
the  first  order  form  of  the  transport  equation.  In  1974, 
Kaper,  Leaf  and  Lindeman  examined  the  application  of  a  finite 
element  projection  of  Lagrange  polynomials  and  surface  har¬ 
monics  to  a  variational  form  of  the  second  order  transport 
equation  in  x-y  geometry.  They  concluded  that  this  approach 
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was  inefficient  and  had  many  disadvantages  in  comparison  to 
the  standard  solution  techniques  of  discrete  ordinates  and 
Monte  Carlo.  Subsequently,  Miller,  et  al.  applied  finite 
elements  to  this  variational  form  of  the  transport  equation 
and  demonstrated  the  mitigation  of  ray  effects.  Pitkaranta 
and  Silvennoinen  also  used  a  phase  space  projection  of 
finite  elements  to  solve  eigenvalue  problems  in  spherical 
geometry;  whereas,  Lillie  and  Robinson  used  a  discrete 
ordinate  finite  element  technique  and  a  reduced  functional 
to  solve  two-dimensional  (x,y)  transport  problems. 

Many  other  researchers  have  used  projection  methods  to 
solve  particle  transport  problems  although  most  of  them 
did  not  explicitly  classify  their  approximate  solution 
techniques  as  projection  methods.  Some  of  those  who  did 
were  Miller  and  Reed  (1975)  who  carried  out  a  discrete 
ordinate  to  spherical  harmonic  conversion  to  mitigate  ray 
effects  and  Reed  (1972)  who  also  developed  a  discrete 
ordinate-spherical  harmonic  solution.  Others  include 
Morel  (1981)  who  used  a  collocation  method  with  Lagrange 
polynomials  to  solve  the  first  order  form  of  the  transport 
equation. 

However  the  most  widely  used  methods  for  solving  parti¬ 
cle  transport  problems  are  still  the  discrete  ordinates  (SR) 
method  along  with  Monte  Carlo  and  the  Pn  method.  These  meth- 
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ods  are  well  documented  in  the  literature  and  have  been 
used  extensively  to  solve  general  problems  (Bell  and 

I  Glasstone,  1971).  There  are  many  S^,  Monte  Carlo  and 

diffusion  production  codes  available  in  the  areas  of 
neutral  and  charge  particle  transport.  Recently  Fletcher 
(1983)  has  numerically  developed  a  high  order  Pn  projection 
method  which  uses  finite  elements  to  solve  the  first  order 
transport  equation.  Nonetheless  and  despite  the  success 
and  robustness  of  these  solution  techniques  there  is  a 
need  for  alternate  and  complementary  solution  methods  for 
solving  particle  transport  problems.  This  is  especially 
true  for  those  problems  where  discrete  ordinates,  Monte 
Carlo  and  diffusion  theory  have  disadvantages  and  limi¬ 
tations  . 

Research  Objective 

The  purpose  of  this  research  is  to  develop  a  method 
for  solving  general  particle  transport  problems.  This  meth¬ 
od  should  serve  as  a  complementary  alternative  to  the 
standard  solution  techniques  of  discrete  ordinates,  Monte 
Carlo  and  the  Pn  method.  The  objectives  of  this  study  are 
grouped  into  three  main  areas .  These  were  to  first  provide 
a  theoretical  and  mathematical  development,  followed  by  a 
numerical,  and  then  a  computational  treatment. 
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In  the  first  area  a  theoretical  and  mathematical  formu¬ 
lation  of  the  problem  is  required.  A  mathematical  formalism 
based  on  a  variational  principle  and  the  positive  definite 
self-adjoint  second  order  forms  of  the  transport  equation 
will  be  provided.  The  boundary  conditions  for  most  physical 
problems  of  interest  are  dealt  with  by  this  classical  minimi¬ 
zation  (variational)  principle.  This  development  which  in¬ 
cludes  general  boundary  conditions ,  anisotropic  sources  and 
scattering  has  properties  which  are  amenable  to,  and  would 
facilitate  standard  numerical  solution  techniques. 

The  numerical  development  requires  an  approach  which  is 
consistent  with  and  takes  advantage  of  the  mathematical  for¬ 
mulation.  A  finite  element  projection  technique  based  on 
spherical  harmonics,  piecewise  polynomials,  and  a  piecewise 
constant  (multigroup)  energy  dependence,  met  this  criteria. 
The  general  numerical  treatment  includes  a  seven  dimensional 
time  and  energy  dependent  problem  with  three  spatial  and  two 
angular  or  direction  variables.  The  treatment  of  boundary 
conditions  is  based  upon  an  equivalent  classical  Raleigh- 
Ritz  solution. 

In  order  to  demonstrate  the  potential  of  this  solution 
approach  in  solving  particle  transport  problems  a  computer 
implementation  of  this  method  is  required.  The  application 
to  a  problem  which  is  of  interest  in  the  area  of  radiation 
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physics  was  chosen.  In  this  computational  treatment  a 
validation  of  the  method  through  a  solution  of  the  air- 
over-ground  problem  is  provided.  Specifically,  this  problem 
is  modeled  in  cylindrical  (r,z)  geometry,  with  an  exponen¬ 
tially  varying  atmosphere,  anisotropic  scattering,  first 
scatter  collision  sources  and  with  the  air-ground  interface 
included.  In  this  demonstration  a  steady- state  one  group 
(monoenergetic)  solution  is  achieved  and  examined.  However, 
a  solution  of  the  energy  and  time  dependent  problem  can  be 
obtained  from  a  straightforward  extension  of  this  steady- 
state  monoenergetic  solution. 

Sequence  of  Presentation 

The  organization  of  this  report  covers  the  three  main 
areas  of  research  by  presenting  a  theoretical,  mathematical 
and  numerical  formulation  of  the  general  problem  and  by  a 
specific  computer  implementation  to  the  air-over-ground 
problem.  In  Chapters  II,  III  and  IV  the  mathematical  model 
is  presented  and  discussed  along  with  an  overview  of  finite 
elements  and  the  second  order  forms  of  the  transport  equa¬ 
tion.  In  Chapter  V  variational  principles,  functionals 
and  boundary  conditions  are  presented  and  discussed.  We 
discuss  projection  methods  and  weak  forms  and  outline  the 
general  numerical  method  and  solution  strategy  in  Chapters 
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VI  and  VII.  An  application  to  the  air -over -ground  problem 
which  includes  a  computer  solution  and  results  is  given 
in  Chapter  VIII.  Finally,  Chapter  IX  contains  conclusions 
and  recommendations  for  future  extensions  of  this  work. 
Here  we  also  discuss  applications  to  other  problems  of 
engineering  and  physical  interest. 


CHAPTER  II 


I 

THE  MATHEMATICAL  MODEL 


The  usual  mathematical  description  of  particle  transport 
is  given  by  the  transport  equation  (Duderstadt  and  Martin, 
1971).  This  equation  is  an  expression  in  a  seven  dimen¬ 
sional  phase  space  of  particle  conservation  in  terms  of 
particle  distribution  functions  (number  densities).  The 
transport  equation  models  general  particle  transport  pro¬ 
cesses.  However,  it  can  be  modified  and  used  to  model 
special  neutral  and  charge  particle  transport  problems. 

These  modificiations  occur  primarily  in  the  collision  terms 
(Schmidt,  1979). 

The  transport  equation  along  with  the  relevant 
boundary  and  interface  conditions  is  mathematically  a 
boundary  value  problem.  This  equation  which  is  the  starting 
point  of  our  solution  strategy  is  outlined  and  examined  in 
the  following  sections.  We  also  briefly  mention  the  main 
solution  methods  for  solving  particle  transport  problems. 
Except  for  the  Monte  Carlo  method  these  solution  techniques 
are  based  primarily  upon  the  standard  and  usual  solution 
approaches  for  solving  initial  and  boundary  value  problems. 
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The  Transport  Equation 

A  derivation  of  the  general  transport  equation,  by 
imposing  particle  conservation  and  balance  on  an  elemental 
volume  in  phase  space  can  be  found  in  Duderstadt  and  Martin 
(1979).  This  model  which  includes  external  macroscopic 
forces  and  long-range  collisions  can  be  written  as 


|^(r,v,t)  +  v-Vn(r,v,t)  +  |(r,t) • ?vn(r,v,t)  =  +^q(r,v,t) 


(1) 


where 

/v 

r 

A 

V 

t 


n(r ,v, t) 


A  /V 

E(r , t) 

£<£, t) 


the  spatial  position  vector 

the  particle  velocity  vector 

the  time  variable  with  respect  to 
some  initial  or  starting  time 

unnormali2ed  particle  phase  space 
density  or  probability  distribution 
function 

spatial  gradient  operator 
velocity  gradient  operator 

^[E(r,t)  +  v  x  B(r,t)J  =  particle  accelera¬ 
tion 

electric  field 
magnetic  field 
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rate  of  change  of  particle  distribution 
as  a  result  of  collisions.  This  tern 
is  usually  called  the  collision  operator 

external  particle  source  and  intrinsic 
particle  sources  due  to  fissions,  (n,2n) 
reactions,  etc. 

Equation  (1)  can  also  be  derived  from  the  BBGKY  heir- 
archy  of  equations  (Krall  and  Trivelpiece,  1973).  The 
BBGKY  hierarchy  is  a  system  of  equations  whereby  the  distri¬ 
bution  function  of  a  many  body  system  is  expressed  in  terms 
of  functionals  of  single  particle  density  distributions. 

A  density  expansion  to  first  order  when  carried  out  produces 
the  transport  equation. 

If  we  consider  a  solution  to  Eq  (1)  at  long  times 
and/or  at  large  distances  from  sources  or  boundaries  the 
problem  can  be  formulated  as  an  initial  value  problem 
(Duderstadt  and  Martin,  1979).  In  such  instances  details 
of  the  initial  and/or  boundary  conditions  can  be  ignored 
and  the  time  relaxation  or  slowing  down  to  an  equilibrium 
or  quasi-equilibrium  state  is  usually  described  and  studied 
as  an  initial  value  problem.  However,  if  we  are  interested 
in  particle  distributions  in  the  vicinity  of  sources  and 
boundaries  we  must  construct  our  solution  in  terms  of  a 
boundary  value  problem.  Here  we  seek  and  develop  a  solution 
to  Eq  (1)  as  a  boundary  value  problem. 


3n\ 
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To  further  develop  the  transport  equation  and  the 
precise  model  which  will  be  used  we  seek  simplifications 
|  to  Eq  (1).  First  we  limit  the  collisional  processes  to 

those  which  are  uncorrelated,  and  instantaneous  or  localized 
in  space.  In  a  sense  we  are  limited  to  collision  events 
i  which  are  well  separated  in  terms  of  mean  free  paths,  i.e., 

the  particle  mean  free  path  must  be  larger  than  the  range 
of  interaction  forces.  Then  we  can  write  the  collision 
operator  as 


*  |  v  j  o,.  (r  ,v)n(r  ,v,  t) 


(2) 


where  d3v  is  an  elemental  volume  in  velocity  space  and 


A  A 

ot(r,v)  =  probability  per  unit  distance  travelled 
that  a  particle  at  position  r  and  with 
velocity  v  will  have  an  interaction,  and 

a  (r,v'-*-v)  =  probability  per  unit  distance  travelled 
°  that  at  position  r  a  particle  with  velocity 

v'  will  undergo  a  collision  and  produce 
a  particle  with  velocity  v;  and 

|v|  =  particle  speed  (a  scalar) 


At  this  point  we  can  rgnore  all  collective  effects  and 
if  electric  and  magnetic  fields  along  with  the  initial  and 
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boundary  conditions  are  given  we  can  now  attempt  to  solve 
Eq  (1)  directly.  However,  this  is  still  a  formidable  task 
and  we  usually  make  further  simplifications.  If  we  cannot 
ignore  collective  effects  and  make  these  simplifications,  as 
is  the  case  in  a  dense  plasma  with  a  small  debye  length, 
then  we  must  attempt  to  solve  Eq  (1)  and  determine  the 
electric  and  magnetic  fields  self-consistently  using  Max¬ 
well's  equations. 

To  obtain  the  usual  neutral  particle  (Boltzmann) 
transport  equation,  we  further  ignore  all  macroscopic 

<t-i 

forces  and  delete  the  force  term  — • Vvn(r ,v,t) .  This  is  a 
reasonable  approximation  if  there  is  no  external  electro¬ 
magnetic  field  or  charge  particles  present,  or  if  the 
external  magnetic  field  is  weak  and  the  particle  distribu¬ 
tion  function  is  isotropic  in  velocity  space.  With  these 
assumptions  we  can  now  write  the  reduced  transport  equation 
(Tran  and  Ligou,  1981)  as 

|^(r,v,t)  +  v-vn(r,v,t)  +j vjat (r ,v)n(r ,v, t) 

=  y*d3v' |v' |og(r,v’->-v)n(r  ,v' ,t)  +  q(r,v,t)  (3) 

The  integral  term  in  Eq  (3)  represents  a  source  of  scattered 
particles  from  all  other  velocities  into  the  velocity  space 
v.  This  term  is  usually  called  the  scattering  k__nel. 
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We  are  now  in  a  position  to  make  a  change  of  variables 
which  is  required  in  the  usual  multigroup  formalism,  and  also 
to  write  Eq  (3)  in  terms  of  the  angular  or  phase  space 
flux.  To  do  this  we  must  first  write  the  angular  or  phase 
space  flux  as 


<j>(r,v,  t)  =  |v|n(r  ,v,t) 


and  the  velocity  vector  as 


^  i  ^  I  A 

v  =  |v|ft 


where  ft  is  the  direction  of  particle  travel  and  the  parti¬ 
cle  speed  |v|  can  be  written  in  terms  of  its  energy  as 


-($■ 


Making  the  change  of  variable  from  v  to  E  and  ft  we  can  now 
write  the  transport  equation  in  terms  of  the  angular  flux 
<|>(r,E,ftft)  as 


1  3  d)  A  A  a  a  A  A  A 

|v|3t <r»E,ft,t)  +  ft-V<j)(r,E,  ft,t)  +  ct(r,E)<j>(r ,E,ft,t) 

=  f"  dE’  fdft'  o  (r,E*  ,6'-*-E,ft)<l>(r  ,E'  ,ft',t)  +  q(r,E,ft,t) 

Jq  -/**7T  S 
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Eq  (7)  is  the  standard  time  dependent  transport  equation, 
which,  with  associated  initial  and  boundary  conditions  is 
primarily  a  description  of  neutral  particle  transport. 
However,  depending  upon  the  treatment  of  the  scattering 
kernel  it  can  also  be  used  to  describe  and  solve  certain 
charged  particle  transport  problems  (Przbylski  and  Ligou, 
1981). 

It  is  important  to  note  that  Eq  (7)  is  an  integral- 

differential  equation  with  seven  independent  variables , 

where  we  have  rewritten  the  angular  or  phase  space  flux 

from  a  function  of  (r,v,t)  =  (x,y,z,v  ,v  ,v  ,t)  to  one 

x  y  z 

where  (r,v,t)  =  (r,E,fi,t)  «  (x,y,z,E,0 ,x,t)  ,  and  r  is 
in  a  rectangular  spatial  geometry.  The  source  term 
q(r,E,0,t)  contains  sources  due  to  fissions  and  other 
absorption  processes  which  depend  on  the  angular  flux. 
Furthermore,  and  although  this  is  not  usually  true,  Eq  (7) 
could  be  non-linear  if  the  scattering  cross-section  (a  ) 
is  expressed  as  a  function  of  the  angular  flux. 

The  Boundary,  Continuity  and  Interface  Conditions 

In  order  to  complete  our  description  of  the  problem 
and  formulate  a  numerical  solution  we  must  first  specify 
initial  and  boundary  conditions.  The  initial  condition  is 
simply  given  by 
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(8) 


4>(r,E,n,te)  “  <j>0(r,E,fi)  for  all  r,  E  and  Q 

where  tc  is  some  initial  time. 

The  boundary  conditions  depend  on  the  problem  of 
interest.  In  this  study  we  will  consider  a  set  of  the 
most  common  boundary  conditions .  These  are : 

a.  Vacuum  boundary.  This  is  a  boundary  condition 
where  particles  are  allowed  to  leave  the  surface  but  are 
not  allowed  to  reenter  it.  This  imposes  the  condition 
that  all  incoming  flux  is  zero  at  the  surface.  In  effect 
this  boundary  condition  says  that  the  flux  is  zero  for  all 
incoming  directions  and  it  is  usually  written  as 

<Hr  ,E,fi,t)  =  0  for  all  fi  such  that  ft*  n  <  0  (9) 

s 

A  ^ 

where  rg  denotes  the  boundary  and  n  is  a  unit  outward 
normal  to  the  surface. 

b.  Incident  source.  In  this  case  particles  are 
allowed  to  enter  and  leave  the  system.  The  entering 
particles  represent  an  incident  source  of  particles  on  the 
boundary.  This  boundary  condition  can  be  written  as 

4> (r  ,E,ft,t)  =  q  (r  ,E,ft,t)  for  all  ft  such  that 
s  s  s 

ft*n  <  0  (10) 

where  q„  is  known. 
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c.  Dirichlet  boundary  condition.  The  angular  flux 
can  be  completely  specified  at  the  boundary  as 

4>(rs,E,6,t)  =  <|>s(rs,E,f2,t)  (11) 

where  4>  (r  ,E,S),t)  is  given  and  could  be  a  known  seven 

5  S 

dimensional  function  specified  at  the  boundary,  or  a  con- 

A  A 

stant  in  some  or  all  of  the  independent  variables  (r,E,f2,t). 

d.  Albedo  or  reflecting  boundary.  For  this  boundary 
condition  we  assume  that  some  of  the  particles  leaving  the 
system  are  reflected  back  into  the  system.  This  condition 
can  be  represented  as  (Duderstadt  and  Hamilton,  1976) 

4>(rg,E,-f2,t)  =  a(rg)4>(rg,E,n,t)  (12) 

A 

where  a(r  )  is  the  local  albedo  or  fraction  of  particles 

9 

A 

which  are  reflected.  If  a(r  )  =  1  ,  then  this  is  the  usual 

s 

condition  of  total  reflection  where  all  the  particles  are 
reflected  back  across  the  surface.  It  is  important  to  note 

A  A 

that  the  albedo  could  be  a  function  of  r  ,ft,E  and  t  (Chilton 

s 

et  al. ,  1984).  However,  in  order  to  simplify  the  notation 

A 

we  have  chosen  to  write  it  as  a  function  of  rg  only. 

The  interface  and  continuity  conditions  are  just  ex¬ 


pressions  of  continuity  of  the  angular  flux  in  space  and 


time  (Bell  and  Glasstone,  1970).  They  can  be  written 


as 

lim  4>(r+  sn,ft,E,t  +  s/v)  =  lim  $(r  -  s£2,6,E,t  -  s/v)  (13) 
s-*-o  s-+o 

where  s  represents  a  distance  of  particle  travel  along  6 
and  r  could  be  any  point  within  the  problem  domain  to 
include  an  interface.  This  condition  requires  the  continuity 
of  particle  flux  in  space  and  time  for  a  given  energy 
and  direction.  Note  that  the  flux  can  be  discontinuous 
in  velocity  (direction  and  energy) ,  and  it  usually  is 
especially  at  vacuum  boundaries  and  at  interfaces  (Bell  and 
Glasstone,  1970). 

Summary 

Now  we  have  a  complete  and  general  mathematical  de¬ 
scription  of  neutral  particle  transport  processes  and  conser¬ 
vation,  as  it  is  represented  by  the  first  order  form  of  the 
time  dependent  transport  equation.  Further  descriptions 
and  simplifications  can  be  made  with  respect  to  the  time, 
energy,  angular  and  space  dependences,  and  also  with 
respect  to  scattering  and  sources.  However,  with  or  without 
these  additional  simplifications,  the  main  solution  techni¬ 
ques  for  solving  the  transport  equation  with  associated 


boundary  and  continuity  conditions  are  the  Pn  method  or 
diffusion  theory,  discrete  ordinates,  and  Monte  Carlo. 

These  methods  are  well  documented  and  a  complete  and 
detailed  description  of  them  can  be  found  in  the  literature. 
See  Bell  and  Glasstone  (1970),  Greenspan,  et  al.  (1968) 
and  Duderstadt  and  Hamilton  (1976).  These  descriptions 
will  not  be  repeated  here.  Instead,  we  will  now  present 
an  alternate  finite  element  method  for  solving  the  second 
order  forms  of  Eq  (7)  with  the  boundary  and  continuity 
conditions  of  Eqs  (9)  through  (13). 
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CHAPTER  III 


THE  FINITE  ELEMENT  METHOD 


The  finite  element  method  is  a  projection  method  where 
the  problem  solution  is  expanded  in  a  finite  sum  of  piece- 
wise  functions.  The  basic  projection  procedure  is  to  divide 
the  problem  space  into  a  finite  number  of  discrete  sub- 
domains  called  elements  and  to  define  the  solution  in  terms 
of  locally  defined  functions  with  unknown  coefficients. 

This  trial  space  of  piecewise  functions,  which  are  locally 
defined  on  individual  elements  of  the  problem  space,  forms  a 
local  basis  where  the  unknown  and  to  be  determined  co¬ 
efficients  are  usually  called  generalized  coordinates  (Desai 
and  Abel,  1972).  Thus  the  fundamental  concept  of  the 
finite  element  method  is  that  the  solution  can  be  approxi¬ 
mated  by  a  local  basis  and  some  projection  or  approximation 
technique . 

Widespread  use  of  the  finite  element  method  began  in 
the  1950's,  with  its  use  in  solving  structural  problems  in 
the  design  of  aircraft.  Since  then  the  method  has  been 
applied  to  problems  in  such  diverse  areas  as  solid  and  fluid 
mechanics,  heat  transfer,  and  particle  transport.  Its 
versatility  and  advantages  in  solving  problems  with 


differing  material  properties,  complex  and  irregular 
geometries,  and  with  a  variety  of  boundary  conditions  have 
contributed  to  the  success  of  this  solution  approach  in 
solving  general  problems  of  physical  and  engineering  in¬ 
terest.  Its  relation  to  a  classical  Raleigh-Ritz  solution 
and  the  method  of  projections,  coupled  with  a  theoretical 
foundation  in  approximation  theory,  has  assured  convergence 
for  many  problems  (Whiteman,  1975).  Also,  the  usual 
practice  of  using  a  local  basis  of  piecewise  polynomials 
and  standard  element  types  and  shapes,  facilitates  the 
numerical  treatment  and  computer  implementation  of  the 
method.  These  factors  coupled  with  the  availability  of 
large  computers  have  made  the  finite  element  method  success¬ 
ful.  We  will  discuss  some  of  these  factors  in  the  following 
sections . 

Discretization  of  the  Problem  Domain 

When  applying  the  finite  element  method  the  problem 
space  must  be  divided  into  subregions  or  elements  (Hinton 
and  Owen,  1979).  The  numerical  analyst  or  engineer  must 
decide  as  to  the  number,  type  and  shape  of  elements  to  be 
used.  He  must  also  decide  upon  the  number  of  nodes  (inter¬ 
polation  points)  within  each  element,  the  type  of  nodal 
variables,  and  the  type  of  interpolation  functions. 


27 


Usually  for  one -dimensional  problems  with  one  independent 
variable  the  element  is  a  straight  line  with  two  nodes  and 
a  linear  Lagrange  polynomial  basis.  However,  if  a  local 
basis  of  higher  order  Lagrange  polynomials  is  used  then 
the  number  of  nodes  must  correspond  to  the  order  of  the 
polynomials . 

In  two  and  three  dimensions  there  is  a  wider  choice  of 
elements  available.  The  most  common  two-dimensional  elements 
are  the  triangle,  rectangle,  and  general  quadrilateral. 

Again,  the  number  and  types  of  nodes  and  nodal  variables 
depend  on  the  element  type  and  the  order  of  the  interpola¬ 
ting  polynomials.  For  problems  in  three  dimensions  with 
axial  symmetry  we  can  construct  axisymmetric  ring  elements. 
These  elements  which  are  defined  by  two  independent 
variables  are  usually  useful  in  cylindrical  geometry.  They 
are  really  based  on  two  dimensional  elements,  which  then 
allow  us  to  construct  axisymmetric  triangles,  rectangles 
or  quadrilateral  ring  elements.  Other  three-dimensional 
elements  are  the  tetrahedron,  parallelepiped  or  right 
prism,  and  the  general  hexahedron  (Tong  and  Rossettos, 

1977). 

In  addition  there  are  a  large  number  and  variety  of 
element  types  and  shapes  other  than  the  most  commonly  used 
ones.  It  is  also  possible  to  construct  elements  with  curved 


sides.  These  elements  are  very  useful  for  modeling  problems 
with  curved  boundaries .  They  are  usually  constructed  in  a 
local  element  coordinate  system  called  natural  coordinates. 
In  these  coordinates  they  are  straight-edged  elements, 
however,  when  they  are  transformed  back  into  the  problem 
or  global  frame  they  produce  elements  with  curved  sides 
(Whiteman,  1976).  Figure  1  shows  some  two  and  three-dimen¬ 
sional  elements. 

Polynomial  Basis  and  Generalized  Coordinates 

The  most  widely  used  and  accepted  trial  functions  are 
polynomials .  These  functions  are  used  to  represent  the  be¬ 
havior  of  the  solution  within  each  element.  They  also  de¬ 
termine  the  number  of  nodes  within  each  element  and  the  type 
of  nodal  variables.  Depending  on  the  approximating  or  shape 
functions,  we  can  have  elements  with  interior  and  exterior 
nodes  (Oden  and  Reddy,  1976).  Exterior  nodes  are  those 
nodes  which  are  on  the  surface  or  boundaries  of  the  element. 
Interior  nodes  are  located  within  the  element  boundaries. 
Furthermore,  interior  nodes  usually  lead  to  a  computational 
procedure  called  condensation  whereby  the  overall  size  of 
the  problem  matrices  and  equations  is  reduced. 

The  unknown  coefficients  of  the  interpolating  poly¬ 
nomials  are  called  generalized  coordinates.  They  are  not 
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local  coordinate 


global  coordinate 


a.  3-node  triangle 


local  coordinate  global  coordinate 

b.  6-node  triangle 


local  coordinate  global  coordinate 

c.  8-node  quadrilateral 


local  coordinate  global  coordinate 

d.  10-node  tetrahedron 


Figure  1.  Isoparametric  Element  Shapes 


associated  with  actual  nodal  values  of  the  solution  and 
therefore  they  do  not  have  any  direct  physical  meaning. 
However,  the  usual  practice  in  the  finite  element  procedure 
is  to  express  these  element  shape  functions  in  terms  of  nodal 
or  point  values  of  the  solution.  These  are  the  unknown 
coefficients  which  are  to  be  determined  from  the  approxi¬ 
mating  procedure  of  the  finite  element  method.  They  now 
represent  the  solution  at  specified  points  (nodes)  within 
the  problem  domain  (Zienkiewicz ,  1971).  This  expression 
of  the  unknown  polynomial  coefficients  in  terms  of  nodal 
values  gives  a  direct  physical  interpretation  to  the  solution 
and  also  insures  inter-element  continuity. 

The  most  common  and  widely  used  polynomial  basis 
is  the  Lagrange  polynomial.  However,  Hermite  and  higher 
order  spline  bases  have  also  been  used.  We  can  derive  the 
Lagrange  polynomials  in  terms  of  the  solution  values  at 
the  nodes;  whereas,  Hermite  polynomials  allow  this  deri¬ 
vation  in  terms  of  the  solution  and  its  first  derivative. 

In  a  Hermite  or  higher  order  spline  basis  the  unknown  co¬ 
efficients  usually  represent  not  only  solution  values  at 
the  nodes  but  also  solution  derivatives  (Prenter,  1975). 

The  interpolating  Lagrange  polynomial  can  be  written 
as  (Atkinson,  1978) 
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I^Cx)  =  n  ^-x--_x— J  .  m  =  0,1. .  .k-l,k+l. .  .n 


th 

where  L^x)  is  a  product  of  n  terms  and  represents  the  k 
piecewise  polynomial  of  degree  n.  At  x  =  x^  the  polynomial 
LjJx)  is  equal  to  one  and  at  x  =  xm  with  m  f  k  it  is  equal 
to  zero.  Therefore  each  polynomial  has  n  zeros. 

For  a  one  dimensional  element  with  n  points  (nodes) 
we  can  usually  approximate  our  solution  4>(x)  by 


n 

4>(x>  «y;w,)  (is; 

i=0 


where  4>  ^  =  (Mx.^)  and  represents  the  solution  values  at  the 
nodes.  L^(x)  then  represents  our  shape  or  trial  functions. 

An  alternate  shape  function  can  be  written  for  this 
element  in  terms  of  a  monomial  basis  and  generalized  coor¬ 
dinates  as  follows 


4>(x) 


(16 


where  cu  represents  the  generalized  coordinates  or  unknown 
coefficients.  A  direct  procedure  for  deriving  the  Lagrange 


basis  Eq  (15)  from  Eq  (16)  can  be  found  in  Heubner,  1975. 
Derivations  and  expressions  of  Hermite  and  general  spline 
bases  can  be  found  in  the  literature  and  will  not  be  re¬ 
peated  here.  Interpolation  polynomials  in  two  and  three 
dimensions  and  for  various  element  types  and  shapes  can 
also  be  found  in  most  books  on  the  finite  element  method. 

It  is  not  our  purpose  to  present  them  here.  However,  we 
will  now  briefly  discuss  the  concept  of  natural  coordinates 
and  parametric  elements. 

Natural  Coordinates 

A  local  coordinate  system  which  depends  on  the  element 
type  and  shape,  and  has  values  which  range  between  one  and 
minus  one  is  called  a  natural  coordinate  system.  These 
coordinates  are  defined  with  respect  to  element  geometry 
and  with  a  linear  variation  between  nodes .  Each  coordinate 
in  a  n-dimensional  system  has  unit  value  at  one  node  and 
zero  value  at  all  other  nodes.  Also,  the  coordinate  func¬ 
tions  are  usually  normalized  to  one,  so  that  a  sum  of  these 
functions  at  any  point  within  the  element  equals  one. 

Natural  coordinates  provides  the  flexibility  of 
constructing  elements  of  different  types,  shapes,  nodes, 
nodal  variables  and  interpolation  functions.  They  are  very 
useful  in  the  development  of  close  form  integration  formulas 
for  evaluating  the  integrals  of  the  problem  equations.  Also 
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they  are  essential  to  the  construction  of  elements  with 
curved  sides. 


The  development  of  natural  coordinates  for  a  large 
number  of  elements  can  be  readily  found  in  the  literature 
(Tong  and  Rossettos,  1977).  In  one,  two  and  three  dimen¬ 
sions  this  development  and  the  concept  of  natural  coordi¬ 
nates  is  usually  synonymous  with  length,  area  and  volume 
coordinates  respectively.  Furthermore,  a  simple  and 
straightforward  prescription  for  defining  a  natural  coordi 
nate  system  is  based  upon  the  use  of  Lagrange  polynomials. 
For  a  three  dimensional  element  we  can  write  the  element 
coordinates  as 


m 


x  = 


z  = 


E 

xiLi(e) 

m 

E 

y^^n) 

i=l 

m 

E 

z-L^B) 

(17 


(18 


(19 


where  e,  p  and  8  range  between  one  and  minus  one,  and  x^, 
y^,  and  z^  represent  nodal  values  of  the  element  in  the 
global  or  problem  coordinates;  m  is  the  number  of  element 
nodes  and  L.  represents  the  usual  Lagrange  interpolant  of 


Eq  (14).  Then  we  can  write  the  element  approximation  or 
shape  functions  as 


m 

4>(e.n,B)  ^.N. (e,n,B)  (20) 

i=l 

where  <j>^  =  <J)^(e,n,B)  and  represents  the  solution  values 
at  the  nodes.  Also 


Ni(e,n,B)  =  Li(e)Li(n)Li(&) 


(21) 


There  are  many  other  ways  of  deriving  a  natural  coor¬ 
dinate  representation  of  the  element  trial  space.  A  similar 
representation  to  Eq  (21)  for  one  and  two  dimensional  ele¬ 
ments  can  also  be  obtained  using  Lagrange  polynomials. 
However,  in  general,  this  approach  is  limited,  especially 
in  the  construction  of  high  order  cubic  and  quintic  ele¬ 
ments  (Huebner,  1975).  For  these  elements,  other  approaches 
are  used,  which  reduce  the  number  of  interior  nodes.  An 
example  of  this  is  the  family  of  high-order  rectangular 
elements  developed  by  Ergatoudis,  et  al.,  (1968).  These 
elements  have  only  exterior  nodes ,  but  they  also  allow 
the  parametric  representation  of  natural  coordinates. 
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For  problems  with  complex  geometries  and  curved 
boundaries  an  accurate  physical  description  can  be  achieved 
in  the  finite  element  representation  by  the  use  of  para¬ 
metric  elements.  These  are  curved-sided  elements  in  the 
global  or  problem  coordinates.  The  idea  is  to  map  or  trans¬ 
form  straight-sided  elements  in  the  natural  coordinates  to 
curved-sided  elements  in  global  coordinates.  Then,  to 
evaluate  the  problem  equations  and  solve  the  problem  in 
this  global  problem  space.  Unlike  a  solution  in  the  natural 
coordinate  space  where  we  can  use  closed  form  integration 
formulas  we  must  now  use  numerical  integration  techniques. 
Numerical  integration  is  usually  required  because  of  the 
complexity  of  the  integrals  in  the  global  problem  space 
(Segerlind,  1976). 

In  transforming  from  the  natural  or  local  coordinates 
to  the  global  coordinates  we  usually  assume  that  this  trans¬ 
formation  (Jacobian)  exists  and  is  unique.  However,  there 
are  prescriptions  and  checks  to  ensure  that  this  is  indeed 
the  case  (Aziz,  1972).  Furthermore,  in  an  isoparametric 
representation  which  satisfies  certain  completeness 
criteria,  existence  and  uniqueness  is  guaranteed  (Huebner, 
1975). 

There  are  three  basic  categories  of  parametric  elements 
depending  upon  the  number  of  element  nodes  and  the  interpola- 


tion  polynomial.  The  underlying  concept  is  that  the  inter¬ 
polating  polynomial  which  determines  and  maps  the  element 
shape  can  be  different  than  those  of  the  solution  or  trial 
space.  For  an  isoparametric  representation  (or  element) 
these  polynomials  are  the  same.  However  for  a  subpara- 
metric  element  its  shape  is  defined  by  a  lower  order  poly¬ 
nomial  than  those  of  the  trial  space  whereas  a  superpara- 
metric  element  is  defined  by  higher  order  polynomials. 

As  an  example,  an  isoparametric  element  could  be  represented 
by  Eqs  (17)  to  (19)  where  the  L^'s  are  the  same  as  those 
of  Eq  (21). 

Application  to  Transport  Problems 

The  application  of  finite  elements  to  transport  pro¬ 
blems  began  in  the  1970's  with  a  finite  element  solution  of 
the  diffusion  equation  (Kaper,  et  al.,  1972).  For  the  one 
speed  diffusion  equation  the  Galerkin  projection  approach 
produces  a  self-adjoint  system  which  is  equivalent  to  a 
classical  Raleigh-Ritz  solution  (Strang  and  Fix,  1973).  Be¬ 
cause  of  this  equivalence  the  finite  element  method  gained 
widespread  acceptance  as  a  viable  approach  to  solving  par¬ 
ticle  diffusion  problems.  Questions  of  existence,  unique¬ 
ness  and  convergence  of  the  solution  could  be  dealt  with 
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in  the  framework  of  a  variational  principle  and  the 
minimization  of  a  functional,  or  in  terms  of  the  related 
energy  methods  of  Miklin  (1966). 

The  finite  element  method  has  also  been  applied  to 
the  transport  equation.  Kaper,  et  al.,  (1974),  Ohnishi 
(1971),  Ukai  (1972),  Reed  (1973),  Fletcher  (1983)  and  others 
have  used  finite  elements  to  solve  the  first  order  forms  of 
the  transport  equation.  Phase-space  finite  elements  in  both 
the  space  and  angle  variables  have  been  used.  Also,  appli¬ 
cations  of  the  discrete  ordinate  method  in  angle  and  the 
finite  element  method  in  space  have  produced  the  Onetran 
(Hill,  1975)  and  Trident  (Seed  et  al.,  1977)  computer 
codes.  Other  finite  element  applications  based  on  the 
second  order  forms  of  the  transport  equation,  have  also 
been  conducted  (Sanchez  and  McCormick,  1982).  Specific 
examples  of  these  applications  will  be  '^esented  in  the 
next  chapter. 
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CHAPTER  IV 


FINITE  ELEMENTS 
AND  SECOND  ORDER  FORMS 


The  finite  element  method  is  usually  a  mathematical 
and  numerical  technique  for  obtaining  an  approximate  solu¬ 
tion  to  a  large  class  of  problems.  Initially  it  was 
developed  and  used  to  solve  problems  in  stress  analysis 
(Heubner,  1975).  Later,  as  the  mathematical  foundation 
of  the  method  was  established  it  gained  widespread  accep¬ 
tance  and  use  in  solving  a  larger  class  of  problems. 

Finite  elements  are  based  on  projection  procedures, 
however ,  they  can  also  be  used  in  the  Raleigh-Ritz  technique 
of  first  recasting  the  problem  in  an  equivalent  variational 
form  and  then  seeking  a  solution  on  the  basis  of  an  energy 
minimization  principle  (Strang  and  Fix,  1973).  Nonetheless, 
the  original  Raleigh-Ritz  procedure  did  not  include  this 
finite  element  approach.  Instead  the  solution  was  expanded 
in  the  form  of  a  linearly  independent  set  of  global  trial 
functions.  Unlike  the  finite  element  method,  these  global 
functions  were  neither  piecewise  polynomials  nor  were  they 
equal  to  zero  on  parts  of  the  problem  domain. 


Using  this  expansion  of  global  trial  functions,  the 
Raleigh-Ritz  procedure  was  to  find  an  extremum  (minimum 
or  maximum)  of  the  variational  principle  (functional) .  If 
this  linear  combination  of  globally  defined  functions  did 
not  extremize  the  functional,  then  the  class  (or  space)  of 
functions  was  expanded  by  the  addition  of  more  functions. 
This  expansion  of  the  trial  space  continues  until  a  linear 
combination  of  functions  which  is  an  extremum  of  the 
functional  is  obtained  (Miklin,  1964). 

In  the  finite  element  approach  the  problem  domain  is 
divided  into  smaller  regions  or  elements ,  and  the  trial 
functions  are  piecewise  polynomials  which  are  zero  on  parts 
of  the  solution  domain  (a  local  basis).  Also,  the  trial 
space  (number  of  trial  functions)  is  expanded  by  using 
more  elements  and  not  by  the  addition  of  a  new  class  of 
functions.  Because  of  these  differences  the  finite  element 
method  is  more  adaptable  towards  a  numerical  (computer) 
solution  than  the  original  Raleigh-Ritz  procedure. 

For  some  problems  the  classical  Raleigh-Ritz  proce¬ 
dure  is  really  a  Galerkin  projection  method.  For  these 
problems  it  can  be  shown  (Finlayson  and  Scriven,  1966)  that 
the  Galerkin  method  is  equivalent  to  Raleigh-Ritz.  An 
identical  set  of  matrix  equations  and  therefore  the  same 
solution  is  achieved  by  either  method.  Therefore,  for  many 


problems  the  use  of  finite  elements  in  a  classical  Raleigh- 
Ritz  solution  can  be  considered  to  be  a  finite  element 
projection  method.  These  concepts  will  be  described  in 
some  detail  in  Chapter  VI  of  this  report. 

Projection  methods  however,  do  not  usually  include 
the  use  of  a  variational  principle  or  the  calculus  of  vari¬ 
ations.  In  some  problems  a  variational  principle  has  not 
been  developed  or  may  not  exist,  and  therefore,  the 
Raleigh-Ritz  approach  cannot  be  used.  Nonetheless,  in  such 
cases,  projection  methods  can  be  used  to  solve  the  problem. 
Therefore  projection  methods  can  be  extended  to  a  wider 
class  of  problems  than  Raleigh-Ritz. 

The  equivalence  between  Raleigh-Ritz  and  the  Galerkin 
method  usually  exists  for  problems  which  are  self-adjoint. 
This  equivalence  is  useful  when  considering  questions  of  ex¬ 
istence,  uniqueness,  convergence,  error  bounds  and  the 
treatment  of  boundary  conditions  (Miklin,  1964).  Therefore, 
a  mathematical  formalism  which  is  self-adjoint  and  equi¬ 
valent  to  a  classical  Raleigh-Ritz  solution  is  desirable. 
Unfortunately  however,  the  first  order  form  of  the  trans¬ 
port  equation  is  non- self -adjoint  (Duderstadt  and  Martin, 
1979).  This  is  due  to  the  scattering  kernel,  and  streaming 
operator  or  gradient  term  of  Eq  (7),  which  are  non-self- 
adjoint  . 
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In  contrast,  the  diffusion  equation  (Duderstadt  and 
Hamilton,  1976)  which  is  similar  to  a  Sturm- Li ouville  pro¬ 
blem  is  self-adjoint  with  an  equivalent  classical  Raleigh- 
Ritz  solution  (Hardin,  1977).  However,  the  diffusion 
equation  does  not  have  the  detailed  angular  information 
which  is  required  in  the  solution  of  most  problems  of 
general  engineering  and  physical  interest.  For  these 
problems  a  mathematical  model  which  is  self-adjoint  and 
retains  the  detailed  angular  information  of  Eq  (7)  is 
required.  This  model  is  the  second  order  form  of  the 
transport  equation  which  in  the  limit  of  diffusion  theory 
reduces  to  the  diffusion  equation. 

The  Even  and  Odd  Parity  Equations 

In  order  to  obtain  a  mathematical  model  which  is 
self-adjoint  and  would  permit  a  rigorous  treatment  of 
general  boundary  conditions,  the  even  and  odd  parity 
forms  of  the  transport  equation  will  be  developed. 

The  starting  point  of  this  development  will  be  a  simpli¬ 
fied  form  of  Eq  (7)  and  the  associated  boundary  conditions 
of  Chapter  II.  These  simplifications  will  result  in  the 
steady  state  and  energy  independent  or  one  group  transport 
equation.  We  will  begin  our  development  with  a  brief 
discussion  of  the  usual  numerical  strategies  for  solving 
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time  and  energy  dependent  problems.  This  is  the  general 
problem,  presented  in  Chapter  II,  which  we  would  like  to 
solve.  Therefore,  our  intent  is  to  explain  that  the 
following  monoenergetic  steady  state  treatment  can  be 
easily  extended  to  time  and  energy  dependent  problems.  A 
complete  and  general  treatment  of  particle  transport 
requires  this  time  and  energy  dependent  solution. 

The  usual  multigroup  treatment  of  the  energy  depen¬ 
dence  in  particle  transport  problems  is  well-established 
(Bell  and  Glasstone,  1970).  This  multigroup  formalism 
involves  approximating  the  energy  variable  by  piecewise 
constant  functions  on  fixed  energy  intervals  called  energy 
groups.  This  produces  a  discontinuous  approximation  in 
energy  and  a  set  of  coupled  one  group  or  monoenergetic 
problems.  The  multigroup  formalism  is  a  straightforward 
extension  of  the  monoenergetic  development  to  include 
multigroup  fluxes,  scattering,  sources  and  a  set  of  coupled 
one  group  problems.  Therefore,  a  one  group  development 
is  directly  applicable  and  can  be  easily  extended  to  an 
energy  dependent  treatment  of  particle  transport. 

In  the  time  dependent  case  there  are  a  number  of 
approaches  available.  These  approaches  include  the  well- 
known  point  kinetics  model  of  reactor  physics,  and  an 
exponential  time  dependence  (Bell  and  Glasstone,  1970).  For 


diffusion  theory,  the  assumption  of  an  exponential  time 
dependence  produces  a  non-linear  problem  and  involves  an 
eigenvalue  search  for  the  largest  time  eigenvalue 
(Wachpress,  1966).  However,  the  usual  discrete  model  is 
to  represent  the  time  derivative  of  Eq  (7)  by  an  explicit 
forward  or  weighted  difference  (Duderstadt  and  Martin, 
1979).  This  finite  difference  approximation  in  contrast  to 
the  multigroup  energy  formalism  retains  the  basic  steady 
state  treatment.  It  involves  a  set  of  implicitly  coupled 
steady  state  (time  independent)  problems  whereas  the 
multigroup  formalism  produces  a  set  of  implicitly  coupled 
one  group  problems.  Therefore,  energy  and  time  dependent 
problems  can  be  solved  by  a  straightforward  extension 
of  our  steady  state  monoenergetic  treatment. 

To  be  precise,  we  can  discretize  the  time  and  energy 
dependence  of  Eq  (7)  and  write  it  as 


.4+1  a1 
♦b  -  Ab 

vbAt* 
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and  then  rewrite  it  as 


(23) 
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where  the  operator, 


L°  =  fi'Vo  + 


Jtb°  "  fd(i'as 

JkTT 


■Q) 


(24) 


and 


L  =  L  + 


vb  Lti 


(25) 


and 


e1  -  qk+1  +  ** 


vbAt*i 


(26: 


Also  At  £  *  where  11  and  b  represent  the  time  and 


energy  indices  respectively,  and  <j)b,  <j,£+1 ,  atb,  ogb^b  and 
are  the  usual  multigroup  defined  quantities  (Bell  and 


Glasstone,  1970).  Therefore,  Eq  (23)  is  a  set  of  one 
group,  steady  state  equations  with  modified  and  redefined 
sources,  fluxes  and  cross-sections.  In  fact,  these  equa- 


tions  represent  the  one  group  steady-state  form  of  Eq  (7) 
with  modified  sources  as  in  Eq  (26)  and  with  a  modified 
total  cross-section  c?t  defined  as 


vbAtil 


(27) 


The  significance  of  Eq  (23)  is  that  we  have  reduced 
the  problem  to  a  set  of  coupled  algebraic  equations  in 
time  and  energy.  Starting  with  some  initial  condition  we 
can  solve  for  the  time  dependence  by  a  simple  marching 
scheme.  At  each  step  we  must  then  solve  implicitly  a 
set  of  coupled  one  group  equations  for  the  group  angular 
flux.  This  is  the  usual  procedure  for  obtaining  a  complete 
seven-dimensional  solution  to  the  transport  equation 
(Dupree,  et  al.,  1971  and  Hill,  1976).  This  fully  implicit 
scheme  is  unconditionally  stable  (Clark  and  Hansen,  1964). 

Since  we  have  reduced  the  problem  to  a  series  of 
steady-state  one  group  problems  we  will  now  drop  the  ex¬ 
plicit  time  and  energy  dependence  and  proceed  with  our 
development.  However,  it  is  important  to  bear  in  mind 
this  obvious  extension  to  general  time  dependent  and 
multigroup  problems.  Following  the  derivation  of  Kaplan 
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and  Davis  (1967),  Wheaton  (1978),  and  Wills  (1981)  we  can 
write  the  monoenergetic  steady-state  form  of  Eqs  (7)  or 

/V  A 

(23)  in  terms  of  the  -fl  vector,  where  represents  the 
direction  opposite  to  This  gives 


-ft*V<J>(r  ,-fi)  +  0t(r)$(r,-S7) 


i 


dn,c_(r,-^*fi,)<{)(r,fi' )  =  Q(r,-n) 

(28) 


where  in  the  case  of  Eq  (23)  the  source  Q  and  the  total 

/v 

cross-section  ot(r)  are  defined  in  accordance  with  Eqs  (26) 
and  (27),  and  we  have  dropped  the  explicit  time  and  energy 
dependent  notation. 

The  even  and  odd  parity  terms  will  now  be  defined 
as 

¥g(r,fl)  =  %{<k£,6)  +  4>(r,-6)}  (29) 

¥u(r,n)  =  %|<j>(r,fi)  -  4> (r , -Q)|  (30) 

Qg(r  ,fi)  =  %{q(r,£)  +  Q(r  ,-fl)}  (31) 
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QU(r,6)  =  %{Q(r,Q)  -  Q(r,-fi)} 


(32) 


o|(r,n*fl')  =  %|as(r,n*f2' )  +  o^r.-ft-ft')}  (33) 

au( r,S‘fi')  =  %(ae(r,fi-fi')  -  a  (r,-fl*fl’)}  (34) 

S  \  S  S  " 


where 


g/\  /v 

(r.fi)  = 
fu(r  ,f2)  = 

QgU,sb  = 
Qu(£,£)  = 
ag( rtfi*n') 

S 


even  parity  flux 
odd  parity  flux 
even  parity  source 
odd  parity  source 

=  even  parity  scattering  cross-section 


11  A  A  A  , 

cru(r  -  ft  ) 

S 


=  odd  parity  scattering  cross-section 


We  can  now  proceed  as  in  Wills  (1981) ,  to  derive 
the  even  and  odd  parity  second  order  forms  of  the 
transport  equation  which  are 


-srvKu(r )n*V'fg(r,$:2)  +  Gg(r)yg(r,fi)  -  Qg(r,fi)  -  STVKu(r)Qu(r ,fl) 


(35) 
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(36) 


We  can  express  the  odd  parity  flux  in  terms  of  the 
even  parity  flux  as 


^(r.ft)  =  Ku (r)  |qu(r ,fi)  -  crv^r ,6)} 


(37) 


and  the  odd  parity  flux  as , 

4'g(r,fi)  =  Kg(r)  |Qg(r,fi)  -  Q-VH^r^)}  (38) 


and  the  operators  Gg,  Gu,  Ku  and  Kg  are  defined  as 
(Wills,  1981), 
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Here  we  have  followed  the  usual  practice  of  expanding  the 
scattering  cross-section  in  spherical  harmonics  (Y^m> ,  and 
the  *  superscript  means  the  complex  conjugate  function. 
Also,  the  even  and  odd  parity  scattering  cross-section 
coefficients  are  defined  (Wills,  1981)  as 

Ic^(r)  for  l  -  even 

(43) 

0  for  H  -  odd 
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U,''. 

Vr> 


for  l  -  odd 


for  l  -  even 


(44) 


where  o^r)  is  the  Legendre  scattering  cross-section  ex¬ 
pansion  coefficients.  These  coefficients  originate  from 
the  usual  expansion  or  fit  of  the  cross-section  data  with 
Legendre  polynomials.  This  expansion  can  be  written  as 

L 

“s'vEttWV  <«> 

£=0 


where  uQ  is  the  angle  of  scatter  (See  Figure  2)  and  it  is 
normally  assumed  that  scattering  depends  only  on  the  angle 
of  scatter  and  not  on  the  azimuthal  angle  or  the  incident 
direction  (Duderstadt  and  Hamilton,  1976).  The  cross- 
section  expansion  of  Eq  (45)  can  now  be  transformed  or 
expressed  in  terms  of  the  problem  coordinates  by  use  of  the 
addition  theorem  (Bell  and  Glasstone,  1970). 

In  order  to  complete  our  development  of  the  second 
order  forms  of  the  transport  equation,  we  will  now  discuss 
the  associated  boundary  conditions.  These  boundary  condi¬ 
tions  which  were  presented  in  Chapter  II  must  now  be  re- 
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cose0  =  fi.fi'  =  cosine  of  the  angle  of  scatter 
yy'  +  V(1  -  y2 ) (1  -  y'2)  cos  $ 
azimuthal  angle 


Figure  2.  Scattering  Angle 


written  in  terms  of  the  even  and  odd  parity  flux.  So  as 
to  reduce  the  algebra  we  will  now  write  them  as 

¥g(rg,fi)  +  Yu(rs,fi)  =  ct(rs)  £*g(rgffi)  -  yu(rg,n)j 

+  bqg(rg,6)  for  frn  <  0  (46) 

where  the  adjoint  boundary  condition  is 

Ys(rg,n)  +  VU(rs,n)  =  a(rg)  ^g(rg,fi)  +  Yu(fs,fi)j 

+  bq„ (r ,Q)  for  Q'n  >  0  (47) 

s  s 

and  o<a(r  )  s  1  »  and  b  is  equal  to  zero  or  one. 

s 

For  the, 

1.  Vacuum  boundary 
a(r  )  =  b  =  0 

S 

and, 

2.  Incident  source 
«<V  =  0  ;  b  =  1 
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3.  Albedo  condition 
o«*  (rg>  <  1  5  b  =  0 

and, 

4.  Mixed  condition 

o  <  o(rs)  <  1  5  b  =  1 

also  for  the, 

5.  Dirichlet  condition 

Vg(rs,S)  +  ¥u(rs,fi)  =  4>s(rs,fi)  (48) 

and  and  ¥u  can  also  be  found  from  Eqs  (29)  and  (30). 
Finally,  the  interface  and  continuity  conditions  require 
that  the  odd  and  even  parity  flux  be  continuous  in  space 
and  time  for  a  given  energy  and  direction. 

We  have  now  presented  the  second  order  forms  of  the 
transport  equation.  In  this  derivation  the  angular  flux 
is  expressed  in  terms  of  even  and  odd  components  on  the 
unit  sphere.  These  components  are  the  even  and  odd  parity 
fluxes  which  are  just  the  sum  and  difference  of  the  for¬ 
ward  and  adjoint  angular  flux  (Bell  and  Glasstone,  1970). 


It  is  imporatnt  to  note  that  unlike  the  first  order  trans¬ 
port  equation,  this  second  order  formalism  is  positive 
definite  and  self-adjoint.  This  is  because  the  operators 
of  Eqs  (35)  and  (36)  are  positive  definite  and  self-adjoint. 
(Kaplan  and  Davis,  1967). 

The  Application  of  Finite  Element  Projections 
to  Second  Order  Forms 

The  application  of  a  finite  element  projection  method 
to  solve  particle  transport  problems  is  not  a  new  concept. 
Kaplan,  in  1961,  proposed  a  projection  technique  for  solving 
three-dimensional  reactor  problems.  This  technique,  which 
he  called  flux  synthesis,  was  based  upon  a  variational 
Galerkin  type  solution  of  the  diffusion  equation.  Later, 
Kaplan,  et  al.,  (1967)  extended  this  idea  to  solve  the 
transport  equation.  They  based  their  solution  on  a  varia¬ 
tional  principle  whose  Euler-Lagrange  equations  are  the 
second  order  forms  of  the  transport  equation.  They  also 
used  special  ellipsoid  trial  functions  and  a  space-angle 
synthesis  approach  to  solve  steady-state  Milne  problems 
with  isotropic  scattering. 

The  concept  of  space-angle  synthesis  is  really  a  pro¬ 
jection  technique  for  modeling  the  angular  dependence  of  the 
transport  equation.  The  angular  flux  is  represented  by  a 


linear  combination  of  known  functions  and  mixing  coeffi¬ 
cients,  which  are  required  to  satisfy  the  problem  equations 
in  an  approximate  sense  (Natelson,  1968).  In  contrast 
discrete  ordinates  or  the  Sn  method  restricts  particles 
to  move  in  certain  fixed  directions  and  allows  only  a 
finite  number  of  degrees  of  freedom  in  angle. 

In  1973,  this  space-angle  synthesis  projection 
technique  was  used  by  Miller,  et  al.  to  solve  a  one  dimen¬ 
sional  neutron  transport  problem.  Later,  in  1973,  they 
extended  this  solution  method  to  two-dimensional  problems 
in  x-y  geometry  with  isotropic  scattering  and  sources. 

A  four- dimensional  tensor  product  of  linear  and  bilinear 
polynomials  on  a  finite  element  grid,  was  used  in  the 
minimization  of  a  functional  for  the  even  parity  transport 
equation.  Their  approach  involved  the  use  of  both  rectan¬ 
gular  and  triangular  elements  and  the  imposition  of  vacuum 
and  reflective  boundary  conditions .  The  vacuum  boundary 
conditions  were  treated  as  natural  conditions  whereas 
reflection  was  an  essential  boundary  condition  to  be  im¬ 
posed  on  the  trial  space.  As  a  result  of  these  calculations 
they  concluded  that  this  finite  element  solution  mitigated 
the  ray  effect  problem  which  is  present  in  discrete 
ordinates  (Lathrop,  1968). 


Other  researchers  including  Kaper,  et  al.,  (1974), 
Briggs,  et  al.  (1975),  Ackroyd  (1979),  and  Blomquist  and 
Lewis  (1980) ,  have  used  finite  elements  to  solve  the 
second  order  transport  equation.  Kaper,  et  al.  (1974), 
applied  finite  elements  based  on  a  variational  formulation 
to  the  two-dimensional  multigroup  transport  equation. 

They  used  high  order  cubic  finite  elements,  surface  har¬ 
monic  tensors ,  and  linear  Lagrange  polynomials  over  the 
angular  domain.  This  was  coupled  with  a  product  of  piece- 
wise  and  global  Lagrange  polynomials  in  space.  Their 
conclusion  was  that  high  order  finite  elements  were  not 
a  viable  alternative  to  discrete  ordinates. 

Subsequently,  Briggs,  et  al.,  in  1975,  applied  finite 
elements  to  the  variational  problem  of  the  second  order 
one  group  two-dimensional  neutron  transport  equation  in 
x-y  geometry.  They  examined  the  potential  of  this  solution 
approach  to  mitigate  ray  effects  and  concluded  that  ray 
effect  mitigation  was  due  to  the  elliptic  operators  of  the 
finite  element  formalism  whereas  the  discrete  ordinate 
equation  was  hyperbolic.  The  elliptic  nature  of  the 
finite  element  equations  allowed  the  coupling  or  averaging 
of  the  particle  streaming  along  different  paths  (directions). 
They  argued  that  this  elliptic  coupling  which  is  also  pre¬ 
sent  in  the  diffusion  equation,  but  absent  in  the  Sn 
method,  is  essential  for  ray  effect  mitigation. 
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CHAPTER  V 


VARIATIONAL  PRINCIPLES 

Variational  methods  have  played  an  important  role 
in  many  areas  of  mathematical  physics.  Many  physical 
problems ,  especially  in  the  areas  of  thermodynamics  and 
classical  and  quantum  mechanics,  are  connected  with  the  cal¬ 
culus  of  variations.  Variational  calculus  and  in  particular 
Hamilton's  principle,  have  played  a  vital  role  in  des¬ 
cribing  the  dynamics  of  particle  motion  (Lanczos,  1949). 
Unlike  Newton's  laws  of  motion,  which  are  based  upon  the 
vector  quantities  of  momentum  and  force,  Hamilton's 
(variational)  principle  provides  an  alternate  mathematical 
description  in  terms  of  energy,  which  is  a  scalar. 

In  this  variational  model  the  basic  concept  is  that 
the  motion  of  a  mechanical  system  can  be  described  by 
finding  the  minimum  (or  maximum)  of  a  definite  integral 
(functional).  This  minimization  or  variational  problem  is 
dealt  with  by  the  calculus  of  variations  (Courant  and 
Hilbert,  1953).  Using  a  variational  calculus  approach, 
it  is  possible  to  show  the  equivalence  between  a  varia¬ 
tional  problem  and  the  usual  description  in  terms  of  a 
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differential  equation.  In  the  case  of  Hamilton’s  principle, 
which  describes  a  mechanical  system,  the  functional  is  de¬ 
fined  in  terms  of  the  kinetic  and  potential  energy  of  the 
system.  However,  upon  taking  the  first  variation  of 
this  functional,  we  can  derive  the  Euler-Lagrange  equation; 
which,  for  this  problem  is  just  Newton's  equation  of 
motion. 

This  alternate  and  equivalent  formulation  of  a  pro¬ 
blem,  in  terms  of  a  differential  equation;  or  as  a  varia¬ 
tional  problem,  leads  directly  to  classical  approximation 
techniques  for  solving  a  large  number  of  problems.  We 
have  already  mentioned  the  classical  Raleigh-Ritz  solution 
for  diffusion  problems,  and  its  equivalence  to  the  Galerkin 
method.  The  energy  methods  of  Miklin  are  also  in  this  cate¬ 
gory  where  there  is  a  direct  relation  between  projection 
methods  and  Raleigh-Ritz.  In  fact,  this  relationship 
can  be  easily  established  for  problems  which  are  self- 
adjoint  (Strang  and  Fix,  1973). 

An  important  consideration  in  problem  solving  besides 
the  usual  analytical  questions  of  existence,  uniqueness, 
convergence  and  stability  is  the  treatment  of  boundary 
conditions.  For  problems  where  a  variational  approach 
is  available  these  issues  can  usually  be  dealt  with  in  a 
simple  and  straightforward  manner.  As  an  example  of  this, 
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we  can  normally  recast  a  linear  self-adjoint  problem  in 
terms  of  a  variational  principle  and  a  quadratic  func¬ 
tional.  We  can  then  consider  these  issues  within  the 
variational  context. 

In  the  following  sections  of  this  chapter  we  will 
briefly  discuss  linear  self-adjoint  variational  problems. 

We  will  examine  the  treatment  of  boundary  conditions ,  the 
classical  Raleigh-Ritz  solution  strategy  and  applications 
to  the  parity  equations.  Our  primary  intent  in  this 
development  is  to  later  establish  an  equivalence  to  the 
Galerkin  method.  In  doing  so,  we  will  then  be  able  to 
treat  the  analytical  issues  and  boundary  conditions  within 
this  variational  framework. 

The  Model  Problem 

It  is  well-known  that  some  boundary  value  problems 
can  be  written  as  a  variational  problem  (functional) ;  whose 
minimum  is  a  function  which  also  satisfies  the  boundary 
value  problem.  Therefore,  instead  of  seeking  a  solution 
to  the  boundary  value  problem  directly,  we  can  seek  to 
find  the  minimum  of  a  functional.  This  is  usually 
accomplished  by  selecting  linear  combinations  of  functions 
on  some  function  space  (admissible  class  of  functions). 

In  practice,  for  linear  problems  with  quadratic  functionals, 
this  minimization  procedure  produces  a  set  of  coupled 
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linear  algebraic  equations  which  are  not  difficult  to  solve. 
However,  for  general  non-linear  problems  this  minimization 
process  can  produce  serious  computational  difficulties. 

A  functional  is  basically  the  integral  of  a  function, 
whose  arguments  or  independent  variables  are  themselves 
functions,  and  which  assigns  a  real  number  to  each  indepen¬ 
dent  variable  (function).  Therefore,  a  functional  is 
an  expression  which  converts  a  function  into  a  scalar.  As 
an  example  of  this,  and  to  make  these  ideas  more  precise, 
we  will  now  present  a  linear  problem  and  an  associated 
quadratic  functional.  Our  goal  is  to  discuss  in  a  simple 
setting  a  linear  boundary  value  problem  and  the  related 
variational  problem.  Later,  we  will  extend  these  concepts 
to  the  transport  problem. 

We  begin  by  writing  the  classical  three-dimensional 
Sturm-Liouville  problem  in  operator  notation  as 


Lu(r)  =  f(r) 


(50) 


where  the  self-adjoint  operator  (Miklin,  1964) 
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(51) 


Lo  =  -V*(p(r)V<>)  +  q(r)° 


and  V  represents  the  usual  gradient  operator.  The  associa¬ 
ted  boundary  condition  can  be  written  as 


a(rg)Vu(rg) *n  +  b(rg)u(rg)  +  c(rg)  =  0 


(52) 


Eq  (52)  is  the  general  mixed  boundary  condition.  Setting 
a,  b  or  c  equal  to  zero,  we  can  obtain  the  Neumann  and 
Dirichlet  homogenous  and  non-homogenous  boundary  conditions. 

The  quadratic  functional  corresponding  to  the  above 
problem,  Eq  (50),  can  be  written  as  (Mitchell  and  Wait, 

1977) 


where  the  first  integral  in  Eq  (53)  is  over  the  entire  pro¬ 
blem  domain  (volume)  and  the  second  is  a  surface  integral 
on  the  problem  boundaries . 

Eq  (53)  represents  a  variational  problem  which  has 
been  studied  extensively.  It  can  be  shown  (Wouk,  1979)  that 
a  linear  problem  of  the  form  of  Eq  (50) ,  has  a  solution 
which  is  the  minimum  of  the  quadratic  functional,  Eq  (53) . 
Furthermore  the  Euler-Lagrange  equation  of  this  functional 
is  the  Sturm-Liouville  equation.  The  usual  mathematical 
issues  of  existence,  uniqueness  and  convergence  are  also 
dealt  with  by  Wouk  and  Miklin.  They  have  shown  that  for 
the  classical  Raleigh-Ritz  minimization  of  Eq  (53)  exis¬ 
tence,  uniqueness  and  convergence  are  assured. 

Therefore,  a  numerical  solution  to  the  Sturm-Liouville 
problem  can  be  obtained  by  solving  the  equivalent  varia¬ 
tional  problem.  This  solution  exists,  and  can  be  readily 
obtained,  because  the  Raleigh-Ritz  minimization  process 
produces  an  operator  which  is  positive  definite  symmetric. 
Based  on  this  positive  definite  operator,  we  can  then 
prove  and  answer  the  questions  of  existence,  uniqueness  and 
convergence.  Furthermore,  and  because  of  the  self-adjoint 
nature  of  Eq  (50),  we  can  show  that  a  direct  solution 
by  the  Galerkin  method  also  produces  a  positive  definite 
operator.  By  either  the  variational  or  Galerkin  approach 
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we  end  up  with  an  identical  algebraic  problem  (Miklin,  1964). 
The  algebraic  problem  is  positive  definite  symmetric,  which 
insures  that  a  numerical  solution  exists . 

These  concepts  for  the  model  Sturm-Liouville  problem 
are  well-known.  They  have  the  usual  mathematical  defini¬ 
tions  and  meanings;  these  will  not  be  repeated  here. 

However,  theorems,  proofs,  definitions  and  a  detailed 
analysis  of  the  Galerkin  method  and  variational  problem 
can  be  found  in  Wouk  (1979).  The  Galerkin  method  is  also 
presented  and  discussed  in  Chapter  VI.  Presently 
we  will  continue  to  introduce  these  concepts  for  the  model 
problem. 

The  Boundary  Conditions.  An  important  consideration 
in  the  formulation  of  a  variational  problem  is  the  treat¬ 
ment  of  boundary  conditions.  It  is  the  usual  practice  to 
include  in  the  functional  representation  of  a  problem  some 
or  all  of  the  boundary  conditions.  For  most  problems  this 
is  a  difficult  undertaking.  There  are  no  standard  or 
straightforward  procedures  for  converting  a  boundary  value 
problem  into  an  equivalent  variational  problem.  However, 
for  linear  self-adjoint  problems  there  are  some  guidelines. 
These  guidelines  are  based  upon  the  relation  between  the 
variational  problem  and  a  projected  Galerkin  solution. 

(See  the  energy  methods  of  Miklin,  1964).  Nonetheless, 


the  questions  of  boundary  conditions  is  often  a  separate 
and  difficult  issue. 

The  usual  approach  is  to  use  experience,  knowledge 
and  insight  to  define  a  functional  which  treats  many  of 
the  boundary  conditions  as  "natural"  conditions.  Natural 
boundary  conditions  are  those  conditions  which  are 
automatically  satisfied  by  a  solution  of  the  variational 
problem.  For  a  minimum  problem  the  combination  of 
functions  that  minimizes  the  functional  will  also  satisfy 
the  natural  boundary  conditions.  These  conditions  are  a 
priori  eliminated  from  consideration  in  the  solution  of 
the  variational  problem.  Instead,  we  must  consider  those 
conditions  which  are  not  satisfied  by  the  minimizing 
function.  These  are  usually  called  essential  or  principal 
boundary  conditions. 

The  essential  boundary  conditions  are  imposed  by 
requiring  that  the  functions  which  minimize  the  functional 
also  explicitly  satisfy  these  conditions.  We  are  only  al¬ 
lowed  to  choose  linear  combinations  of  functions  in  the 
minimizing  sequence  that  satisfy  all  the  essential  boundary 
conditions.  However,  we  still  must  determine  which  condi¬ 
tions  are  natural  or  essential.  Although  there  are  some 
guidelines  in  this  regard  (Miklin,  1964),  the  usual 
practice  is  to  examine  the  minimum  or  stationary  point  of 
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the  functional.  This  involves  the  calculus  of  variations 
and  the  concept  of  a  variation  which  we  will  now  introduce 
Let  us  suppose  that  u(r)  is  the  function  which  gives 
Eq  (53)  its  smallest  value,  and  u*(r)  is  another  function 
which  is  almost  identical  to  u(r).  On  the  problem  domain, 
we  regard  u*(r)  to  be  infinitesimally  different  from  u(r) . 
Then  we  define  the  first  variation  of  u(r)  to  be 

3u  =  u*(r)  -  u(r)  =  e<Kr)  (5 

and  the  second  variation  to  be 


32u  =  9(3u)  (5 


where  e  tends  to  zero,  and  $(?)  is  an  arbitrary  continuous 
differentiable  function.  Therefore,  the  variation  of  u(r) 
is  a  measure  of  a  very  small  change  in  u(r)  for  a  given 
value  of  r.  In  taking  the  variation  we  consider  the  inde¬ 
pendent  variable  r  to  be  a  constant  whereby 


(56) 


"T 


3r  =  0 


We  can  also  show  (Huebner,  1975)  that 


V(3u(r))  =  3(Vu(r)) 


(57) 


and 


3 /*  u(r)dr  =  f  3u(r)dr 
JR  JR 


(58) 


where  R  represents  some  volume  in  phase  space. 

To  find  the  stationary  point  of  a  functional  I(u), 
we  require  that  the  first  variation  vanishes.  If  the  second 
variation  is  positive  then  the  stationary  point  is  a  minimum: 
if  it  is  a  negative  then  the  stationary  point  is  a  maximum. 
The  second  variational  derivative  therefore  determines 
whether  we  have  a  maximum  or  minimum  variational  problem. 

In  either  case  we  need  to  take  the  first  variation  of  the 
functional  in  order  to  determine  the  conditions  that  u(r) 
must  satisfy  so  as  to  insure  that  I(u)  has  its  extremum. 
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These  conditions  are  just  the  Euler -Lagrange  equation  and 
the  natural  boundary  conditions  for  the  variational  pro¬ 
blem  (Gelfand  and  Fomin,  1963).  The  Euler-Lagrange  equa¬ 
tion  is  the  partial  differential  equation  of  the  boundary 
value  problem.  For  the  functional  of  Eq  (53)  the  Euler- 
Lagrange  equation  and  natural  boundary  conditions  are 
Eqs  (50)  and  (52). 

In  Appendix  A  we  have  derived  general  expressions  for 

the  first  and  second  variation  of  an  arbitrary  first  order 

functional.  In  Appendix  B  we  use  these  expressions  to  show 

that  if  a(r  )  =  p(r  )  ,  Eqs  (50)  and  (52)  are  indeed  the 
s  s 

Euler-Lagrange  equation  and  natural  boundary  condition  of 
Eq  (53).  Furthermore,  it  can  be  shown  that  the  usual  Neu¬ 
mann  boundary  condition  is  also  a  natural  boundary  condi¬ 
tion.  However,  the  Dirichlet  condition  is  essential  and 
must  be  satisfied  by  the  trial  space. 

The  Raleigh-Ritz  Method.  The  usual  numerical  proce¬ 
dure  for  finding  the  function  u(r)  which  makes  the  func¬ 
tional  I(u)  stationary  is  the  classical  Raleigh-Ritz  method 
(Miklin,  1964).  For  the  model  Sturm-Liouville  problem  this 
stationary  point  is  a  minimum.  Therefore,  the  solution 
to  this  problem  requires  that  we  find  the  minimum  of  Eq  (50). 
This  minimum  is  approximated  by  the  Raleigh-Ritz  method. 
Because  of  the  nature  of  the  variational  problem  (Eq  (53)), 


the  Raleigh-Ritz  method  will  provide  us  with  an  upper 
bound  to  the  minimum  of  I(u);  or  at  best,  we  will  compute 
this  minimum  exactly.  However,  in  practice,  this  is 
never  the  case  and  our  approximate  solution  bounds  the 
true  minimum  from  above.  For  some  problems  (in  elasti¬ 
city  theory)  the  minimum  of  the  functional  represents  the 
minimum  potential  energy  of  the  system  and  therefore  it 
is  physically  meaningful.  However,  for  most  problems 
a  physically  pleasing  interpretation  of  I(u)  is  not 
available . 

In  the  Raleigh-Ritz  procedure  we  minimize  I(u) 
with  a  linear  combination  of  functions  called  a  minimizing 
sequence  (Courant  and  Hilbert,  1953).  These  are  really 
trial  functions  with  unknown  coefficients  and  for  the 
finite  element  method,  they  can  be  expressed  in  the  form  of 
Eq  (20) ,  i. e. 


u(f) 


i  =n 


♦^(r) 


(59) 


Then  we  select  the  coefficients  4^  so  that  I(u)  is  a  mini¬ 
mum  by  requiring  that 
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Eq  (60)  produces  a  set  of  linear  algebraic  equations  in 
the  unknown  coefficients  <f>^.  We  can  then  solve  this 
linear  system  and  thus  determine  the  minimizing  function 
of  I(u). 

For  the  model  problem,  Eqs  (59)  and  (60)  when 
inserted  into  Eq  (53)  gives 


2 p (r ) Vu • VNi  +  [q(r)u  -  f(r)]Ni|dr 


+  2<f[b(r  )u  +  c(r  )]N.ds  =  0  osisM  ( 

•/c  ®  SI 


where  we  can  now  solve  for  the  Furthermore,  the 

matrix  problem  of  Eq  (61)  is  positive  definite  symmetric 
(Huebner,  1975)  and  therefore  a  numerical  solution  can 
be  obtained. 


[i] 


The  parity  equations  are  similar  to  our  model  Sturm- 
Liouville  problem  in  that  the  operators  are  positive  de¬ 
finite  and  self-adjoint  (Kaplan  and  Davis,  1967).  Because 
of  this  we  would  like  to  duplicate  our  variational  treat¬ 
ment  of  the  previous  section  for  the  second  order  trans¬ 
port  problem.  Specifically,  we  will  provide  a  variational 
formalism  for  the  even  parity  transport  equation  to  include 
the  general  boundary  conditions  of  Eqs  (46)  to  (48).  We 
have  chosen  to  use  the  even  parity  equation  for  reasons 
to  be  discussed  later.  However,  the  odd  parity  develop¬ 
ment,  except  for  a  change  of  notation,  is  similar  to  that 
for  the  even  parity  equation. 

Variational  principles  for  the  parity  equations 
already  exist  (Kaplan,  et  al.,  1967).  These  principles 
usually  treat  the  vacuum  boundary  condition  as  a  natural 
condition.  However,  a  treatment  and  discussion  of 
general  boundary  conditions  have  not  been  accomplished. 

To  date,  these  variational  principles  have  been  used  to 
primarily  provide  a  theoretical  basis  for  the  Marshack 
boundary  conditions  of  the  Pn  method  (Davis,  1966);  and 
to  develop  numerical  schemes  for  eliminating  ray  effects 
(Briggs,  et  al.,  1975).  The  variational  principle  of 
Vladimirov  (1963)  is  also  in  this  category.  Furthermore  the 


Roussopoulos  principles  of  Cobb  (1971)  and  Selegnut 
(1959)  are  based  on  the  non-self-adjoint  first  order  trans¬ 
port  equation.  These  and  other  variational  principles 
have  also  been  used  to  accurately  compute  integral  or 
average  quantities  such  as  reaction  rates,  cross-sections 
and  eigenvalues  (Duderstadt  and  Martin,  1979). 

We  begin  by  using  experience,  insight,  trial  and 
error,  and  finally  good  judgment  to  write  a  functional  for 
the  even  parity  transport  equation  as  follows 


I(^g) 


V4»g,KU(n-Vg)> 


+  <H,g ,  Gg'i,S> 


-  2<fi* V'Fg,KuOu>  -  2<'i,g,Qg>]dr 


+ 


[<|a-n|Qeb,H'g(rs,^)>]ds 


(62) 


where  <•,*>  represents  the  usual  inner  product  in  angle 
( 4 7T  steradians),  and  |*|  means  the  absolute  value,  and  n 


is  the  outward  unit  normal  to  the  problem  boundaries.  The 
parity  operators  and  sources  have  the  usual  meaning  and 


Qeb  =  a*8(rg,fi>  -  bqg (rg  ,fi)  (63) 


Therefore,  Qe^  is  a  function  of  ,  r,  and  and  we  let 


l-a(fg) 
~  a(fg)  +1 


(64) 


where  o^a(r  )<  1  ,  and  b  and  a  depend  on  the  boundary 
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condition.  They  are  defined  as 

1 .  Vacuum  boundary 

a(rg)  =0  ;  b  =  0  (64a) 

2.  Incident  source 
cc(rg)  =  0  ,•  b  =  2 


(64b) 


3.  Albedo  condition 

o<  a(r  )  i  1  ;  b  =  0  (64c) 

4.  Mixed  condition 

os  a<*8>  s  1  ;  b  =  ^ra(i  )  (64d) 

s 

and 

5 .  Dirichlet  condition 

2 

o<  o(£s)  <  1  5  b  =  o  or  j  (64e) 

s 

where  we  require  that  ¥®(r,D)  be  a  constant 
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in  some  or  all  of  the  independent  variables , 

and  also  satisfy  one  of  the  above  mentioned 

conditions . 

We  now  prove  that  our  postulate  is  correct  and  that 
the  functional  of  Eq  (62)  is  indeed  the  variational  problem 
corresponding  to  the  second  order  (even  parity)  transport 
problem.  This  proof  which  is  an  extension  of  that  for  the 
model  Sturm-Liouville  problem,  is  carried  out  in  Appendix 
C.  There  we  also  show  that  the  boundary  conditions  (1)  to 
(4)  are  natural  conditions  and  that  the  Dirichlet  condition 
is  an  essential  condition. 

It  is  important  to  note  that  with  vacuum  boundary 
conditions  the  functional  of  Eq  (62)  reduces  to  that  of 


Kaplan  and  Davis  (1967)  with  Mar shack  boundary  conditions 
(Davis,  1966).  However,  we  have  now  included  the  albedo, 
incident  source,  mixed  and  interface  conditions  as  natural 
boundary  conditions.  We  have  also  rigorously  shown  that 
this  concept  of  natural  and  essential  boundary  conditions 
is  indeed  valid.  Furthermore,  we  have  developed  a  quadra¬ 
tic  functional  which  can  be  used  to  provide  a  numerical 
solution  to  general  particle  transport  problems.  This 
functional  has  mathematical  and  numerical  properties  which 
are  similar  to  that  of  our  model  problem. 

Therefore,  we  could  seek  a  solution  (minimum)  based 
on  the  usual  Raleigh-Ritz  procedure  and  expect  that  our 
numerical  problem  would  be  positive  definite  and 
symmetric  (Wills,  1981).  Applying  the  Raleigh-Ritz  method 
to  Eq  (62)  we  get  with 

f8(r,fi)  =  ^  ♦^(r.fi)  (65) 

i=0 


[<n*V1'g,KU(Q-VNi)>  +  <4'g,GgNi> 


<f$-VNi,KuQu>  -<Ni,Qg>]dr 


+  ps  [<  |fi*n|  [2aTg  -  bqg]  ,Ni(rg  ,^)>]ds  =  0  (66) 


CHAPTER  VI 


PROJECTION  METHODS  AND  WEAK  FORMS 


Projection  methods  can  be  used  to  provide  approximate 
solutions  to  a  large  number  of  problems.  They  have  been 
used  extensively  to  solve  boundary  value  and  other  problems 
which  are  described  by  differential  and  integral  equations. 
These  methods  are  important  in  that  they  facilitate  the 
recasting  of  general  problems  into  a  form  more  suitable 
for  a  numerical  solution.  Most  projection  methods 
involve  the  use  of  a  contraction  mapping  (Stakgold,  1979). 
The  original  problem  is  usually  projected  unto  a  reduced 
finite  dimensional  subspace  and  a  solution  is  then  sought 
in  this  new  and  restricted  space. 

Projection  methods  include  such  techniques  as 
collocation,  Galerkin,  the  method  of  moments  and  least 
squares  (Golberg,  1978).  These  methods  are  sometimes 
called  the  method  of  weighted  residuals  (Finlayson  and 
Scriven,  1966).  However,  they  all  begin  with  the  use  of 
a  projection  operator  defined  as  (Atkinson,  1976) 


Pn(x)x  =  x  for  all  x  e  Xn 


(67) 


where  Pn  is  a  bounded  projection  operator  from  X  unto  Xn. 
X  is  a  Banach  space  and  Xn  a  finite  dimensional  subspace 
of  X.  Then 


(68) 


and 


IIPJI  =  I  1**11  I  I pn I  i 2  (69) 

where  11*11  represents  the  usual  definition  of  a  norm  and 
Eq  (69)  implies  that  |  JPn  |  |>1  (Atkinson,  1976). 

Therefore,  a  projection  operator  maps  every  element 
of  a  general  space  X  into  Xn,  a  subspace  of  X.  It  also 
maps  each  element  of  Xnunto  itself.  This  projected  space 
Xn  is  usually  the  n-dimensional  space  in  which  we  seek  to 
find  a  solution. 

If  P_  is  a  projection  operator  then  (I  -  P  )  is  also 
11  n 

a  projection  operator  where  I  represents  the  identity  opera¬ 
tor  (Wouk,  1979).  Furthermore,  in  a  Hilbert  space  we  can 

1 

■  1 
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introduce  the  concept  of  orthogonal  projections  where  the 
projection  operator  involves  the  inner  product  of  two 
functions.  An  orthogonal  projection  operator  Pn  is 
symmetric  and  orthogonal  to  (I  -  P  ) .  Therefore, 

<Pnx, (I  -  Pn)y>  =  o  for  x,yeX  (70) 


where  (I  -  Pn>  is  a  projection  operator  unto  the  orthogonal 
complement  of  XR.  Finally,  for  an  orthogonal  projection  Pn 


ll*ll2  -  l|rnx||!  +  1 1  <i  -  VXH“  <71) 


and 


I 


=  1 


i 


(72) 


The  projection  method  for  solving  a  problem  of  the 


form 


Ly(x)  =  f  (x ) 


(73) 


is  to  approximate  y (x)  on  a  finite  dimensional  space  by 


(74) 


and  define 


R(x)  =  Lyn(x)  -  f (x) 


(75) 


where  R(x)  would  be  identically  equal  to  zero  if  y(x)  = 
Yn(x).  However,  since  in  general  this  is  not  true,  we 
try  to  choose  yn(x)  so  that  R(x)  is  small.  We  accomplish 
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this  by  making  the  projection  of  R(x)  unto  Y^(x)  equal 
to  zero,  i.e. , 


PnLyn(x)  "  Pnf(x)  =  PnR(x)  =  0  (76) 


where  yn  e  Yn  and  Yn  is  the  finite  dimensional  subspace 
of  the  Banach  space  Y.  Therefore,  the  projection  method 
for  solving  Eq  (73)  is  to  pick  yn(x)eYn(x)  in  accordance 
with  Eq  (76)  such  that  the  components  in  yn(x)  of  Lyn(x) 
agree  with  f(x). 

Eq  (76)  represents  a  set  of  coupled  algebraic  equa¬ 
tions.  However,  before  we  can  seek  a  numerical  solution 
we  must  first  choose  a  projection  operator  P  .  There  are 
a  number  of  possibilities,  but  the  usual  selections  lead 
to  the  well-known  methods  of  collocation,  Galerkin,  least 
squares  and  the  method  of  moments.  In  the  collocation 
method  we  require  that  the  residual  vanish  at  selected 
node  points.  Thus  we  determine  the  otj's  of  Eq  (74)  by 
requiring  that 


R(xi)  =  0  l*iSN 


(77) 


* 
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or  equivalently  by 


f 


R(x)w^(x)dx  =  0 


IsisN 


(78) 


where  w^(x)  =  <S(x  -  x^) .  The  Dirac  delta  function  w^(x) 
represents  the  test  or  weight  function  of  the  collocation 
method  (Finlayson  and  Scriven,  1966).  Therefore,  we  have 
defined  or  chosen  our  projection  operator  with  respect  to 
the  weight  function  w^x)  and  in  accordance  with  Eq  (78)  or 
(77). 

For  the  methods  of  Galerkin,  least  squares  and  the 
method  of  moments  we  can  also  require,  as  in  the  case  of 
collocation,  that  the  residual  be  orthogonal  to  some  weight 
function.  This  in  turn  leads  to  the  algebraic  problem 
(Eq  (78)),  and  with  a  proper  choice  of  the  subspace  Yn, 
to  a  solution  for  the  a's  of  Eq  (74).  The  weight  functions 
are  selected  as  follows  (Prenter,  1975). 

a.  Galerkin.  The  residual  is  required  to  be 
orthogonal  to  the  trial  space.  Here  we  take  the  weight 
functions  to  be  the  same  as  the  trial  functions  of  Eq  (74), 
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w^x)  =  4)^  <x) 


b.  Least  Squares.  We  minimize  the  two-norm 
of  the  residual  by  choosing  our  weight  functions  to  be 


wi(x) 


9R(x) 

9a^ 


(80) 


c.  Method  of  Moments.  The  residual  is  made  ortho¬ 
gonal  to  the  monomial  basis  |xn| .  We  choose  the  weight 
functions  to  be 


wi(x) 


x1  oSi*N 


(81a) 


Weak  Forms 

The  most  widely  used  projection  methods  are  colloca¬ 
tion  and  Galerkin.  However,  the  Galerkin  method  is  a 
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generalization  of  the  Raleigh-Ritz  method  for  problems 
which  are  positive  definite  and  self-adjoint.  For  these 
problems  we  can  obtain  a  Raleigh-Ritz  solution  without 
having  to  develop  a  variational  formalism.  Furthermore, 
the  Galerkin  method  is  more  general  than  this  in  that  it 
can  be  applied  to  problems  which  are  not  positive  definite. 

In  applying  the  Galerkin  method  to  two  point  boundary 
value  problems  the  normal  procedure  is  to  use  integration 
by  parts  to  produce  the  Galerkin  weak  form  (Strang  and  Fix, 
1973).  This  procedure  lowers  the  continuity  restriction 
of  the  trial  space  and  allows  for  the  inclusion  of  natural 
boundary  conditions.  We  can  then  establish  a  precise 
equivalence  with  the  Raleigh-Ritz  variational  problem. 

This  equivalence  is  important  for  reasons  which  we  have 
already  mentioned. 

For  the  Sturm-Liouville  boundary  value  problem  of 
Eq  (50)  we  can  use  the  divergence  theorem  to  derive  the 
weak  form.  Assuming  a  trial  space  in  the  form  of  Eq  (59), 
the  Galerkin  projection  can  be  written  as 


j^t-7- p(f)V-u(f)  +  q(t)  -  f(f)]N.(f)df  =  0 

for  o*i*M 


(81b) 


Then,  applying  the  Gauss  divergence  theorem  to  the  first 
term  of  Eq  (81b) ,  we  get 


+  [q(r)u  -  f(r)3N^jdr  -  ^Ep(rg)Vu(rs)*n]N^ds  =  0  l£i£M 

(82) 


and  using  Eq  (52)  we  can  rewrite  Eq  (82)  as 


^|p(r)Vu-VNi  +  [q(r)u  -  f^JN^dr  +  ^(M£s)u(rg)  +  c(^g)]Nids  -  0 

l<i*M  (83) 


Eq  (83)  is  identical  to  Eq  (61).  We  have  therefore  estab¬ 
lished  the  equivalence  of  the  Raleigh-Ritz  solution  and  the 
Galerkin  weak  form.  The  boundary  conditions  of  Eq  (52)  can 
now  be  treated  as  natural  boundary  conditions.  Furthermore 
all  boundary  conditions  for  a  Galerkin  weak  form  solution  to 
Eq  (50)  can  be  treated  within  the  context  of  the  variational 
problem,  as  either  natural  or  essential  conditions. 


CHAPTER  VII 


THE  FINITE  ELEMENT  PROJECTION  METHOD 


In  Chapter  V  we  presented  a  variational  formulation 
of  the  second  order  transport  problem.  However,  we  would 
like  to  provide  a  numerical  treatment  that  is  straight¬ 
forward  and  avoids  the  explicit  use  of  a  variational  princi¬ 
ple.  The  projection  methods  which  we  discussed  in  Chapter 
VI  can  be  used  to  provide  this  alternate  numerical 
treatment.  Furthermore,  for  the  Galerkin  method  we  can 
establish  an  equivalence  to  the  variational  problem.  The 
purpose  of  this  Chapter  is  to  establish  this  equivalence  and 
to  outline  a  finite  element  Galerkin  method  for  solving 
general  particle  transport  problems. 

Our  intent  is  to  provide  a  numerical  solution  to  the 
even  parity  second  order  transport  equation.  In  doing  so 
we  will  treat  the  boundary  conditions  within  the  context 
of  a  classical  Raleigh-Ritz  solution,  as  essential  or 
natural  conditions.  This  classical  solution  was  presented 
in  Chapter  V,  where  we  discussed  in  detail  the  concept 
of  natural  and  essential  boundary  conditions. 

We  begin  by  referring  to  Chapter  IV  where  the  angular 
flux  was  written  in  terms  of  odd  and  even  components 
as 
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<Kr,fi)  =  V%(r,n)  +  fu(r,fi) 


(84) 


where  the  even  and  odd  parity  flux  and  'f'11  can  be  obtained 
by  solving  Eqs  (35)  and  (36).  Therefore,  a  solution  to 
these  equations  provide  a  solution  for  the  particle  angular 
flux.  Furthermore,  because  and  are  even  and  odd 
on  the  unit  sphere  they  represent  the  scalar  or  reaction 
rate  flux  and  the  particle  current  respectively  (Miller 
et  al.,  1973).  For  most  problems  we  are  interested  in  an 
accurate  description  of  the  scalar  flux  and  thus  a 
solution  to  Eq  (35)  is  all  that  is  required.  Because  of 
this ,  our  solution  will  be  presented  in  terms  of  the  even 
parity  transport  equation.  Nonetheless  the  solution  stra¬ 
tegy  can  be  applied  to  either  problem. 

The  Galerkin  Weak  Form 

Our  solution  is  based  on  a  finite  element  Galerkin 
solution  of  the  even  parity  transport  equation.  We  have 
chosen  this  approach  because  solutions  to  Eq  (35)  by  use 
of  a  variational  principle  or  the  Galerkin  method  are 
equivalent  (Wills,  1981).  To  show  this,  we  can  expand  the 
even  parity  flux  as  in  Eq  (65).  Using  the  Galerkin 
projection  operator  we  can  then  write  Eq  (35)  as 


^[<-6'V KUfi-V't'g,Ni>  +  <Gg4'g,Ni>  -  <Qg,Ni> 

+  <fi* VKuQu,Ni>]dr  =  0  (85) 

We  can  now  use  the  divergence  theorem  (Kaplan  and  Davis, 
1967) 


^<f2-Vf,g>dr  -  -^<f,Q-Vg>dr  +  j£<  (6 -h)f  ,g>ds 


(86) 


to  include  the  boundary  terms  of  Eqs  (46)  and  (47).  This 
then  allows  us  to  rewrite  Eq  (85)  as 


2^[<f2*VH'g,KU(f2'VNi)>  +  <fg,GgNi>  -  <fl*VNi,KuQu> 

-  <N.,Qg>]df  +  j£<|6*n|  [2a4'g  -  bqg]  ,N.(fs  ,^)>dS  =  0  (87) 

where  a  and  b  are  defined  by  Eqs  (64a)  to  Eq  (64e). 

Eq  (87)  is  a  set  of  coupled  algebraic  equations  which 
represent  the  Galerkin  weak  form  of  the  even  parity  trans¬ 
port  equation.  It  is  identical  to  the  Raleigh-Ritz  system 
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of  equations,  Eq  (66).  Therefore,  we  have  established  the 
equivalence  of  the  Raleigh-Ritz  solution  and  the  Galerkin 
weak  form.  Furthermore,  we  can  now  consider  the  boundary 
conditions  of  Eqs  (46)  and  (47)  to  be  natural  conditions 
and  the  Dirichlet  condition  to  be  essential. 

This  treatment  of  the  boundary  conditions  and  a 
numerical  solution  of  Eq  (87)  is  our  finite  element  pro¬ 
jection  method.  This  zig-zag  between  a  variational  prin¬ 
ciple  and  the  Galerkin  method  allows  us  to  rigorously 
treat  the  boundary  and  interface  conditions  as  natural  or 
essential.  In  effect,  we  need  a  variational  principle 
in  order  to  explain  our  treatment  of  the  boundary  terms. 

We  also  need  a  variational  principle  in  order  to  determine 
which  boundary  conditions  are  essential  or  natural. 


Solution  Strategy 

At  this  point  of  the  development  we  have  reduced  the 
problem  to  a  set  of  algebraic  equations  once  we  define  the 
trial  basis,  Eq  (65).  The  problem  as  represented  by  Eq  (87) 
is  therefore  a  matrix  problem  for  the  flux  expansion  co¬ 
efficients  (c()i) .  Furthermore,  for  a  general  time  dependent 
multigroup  treatment  (see  Chapter  IV)  we  have  to  solve  a 
number  of  these  matrix  problems,  which  are  then  coupled 
in  time  and  energy.  However,  a  solution  to  Eq  (87)  is 
still  a  formidable  task. 


Nonetheless,  we  can  develop  a  solution  strategy  which 
is  based  on  the  special  properties  of  the  problem  matrix. 
These  properties  together  with  the  parity  (odd,  even) 
characteristics  of  the  even  parity  flux  should  ease  some  of 
the  computational  difficulties.  Also,  by  a  judicious  choice 
of  the  trial  functions  we  can  effectively  treat  the  spatial 
and  angular  dependences . 

Trial  Functions.  For  the  most  general  particle  trans¬ 
port  problems  the  steady  state  one  group  problem  of  Eq  (87) 
is  a  five  dimensional  problem.  The  particle  flux  has 
three  independent  spatial  variables  in  r  and  two  angular 
variables  in  fi.  A  spatial  solution  is  usually  carried  out 
in  rectangular,  cylindrical  or  spherical  geometries.  The 
gradient  or  streaming  operator  (fi’V'O  for  these  orthogonal 
geometries  are  well-known  and  can  be  found  in  Bell  and 
Glasstone  (1970).  Our  intent  is  to  discuss  a  general 
solution  strategy  applicable  to  these  geometries.  Here  we 
will  omit  some  of  the  specific  algebraic  details.  However, 
we  will  address  some  of  them  in  cylindrical  coordinates  and 
as  they  apply  to  the  air-over -ground  problem  which  is  dis¬ 
cussed  in  Chapter  VIII. 

In  Chapter  II  we  wrote  the  particle  velocity  or 
vector  dependence  in  terms  of  particle  energy  and  the 


A 

direction  vector  fl.  This  vector  is  related  to  the  spatial 
problem  coordinates  and  is  usually  specified  by  two 
independent  angle  variables  y  and  x  (Duderstadt  and  Martin, 
1979) .  These  angles  are  defined  with  respect  to  the 
particular  problem  geometry.  In  rectangular  geometry  they 
are  defined  as  in  Figure  3  where  u  =  cos 6  and  8  is  the 
angle  formed  by  ft  and  the  z-axis.  x  is  the  angle  between 
the  x-z  plane  and  the  plane  formed  by  the  ft  vector  and  the 
z-axis . 

By  expressing  the  problem  in  an  orthogonal  spatial 
coordinate  as  in  Figure  3  we  can  easily  derive  expressions 
for  the  streaming  operator.  We  can  also  derive  expressions 
for  the  surface  normal  term  fi-h  (Wills,  1981).  Furthermore, 
we  can  easily  demonstrate  the  odd-even  properties  of  our 
second  order  formulation.  These  properties  are  important 
in  determining  the  scalar  flux  and  particle  current.  They 
also  show  that  this  second  order  formulation  is  just  a 
separation  of  the  odd  and  even  components  of  the  angular 
flux.  The  even  parity  equation  and  flux  carries  the 
even  information  which  is  important  in  determining  the 
scalar  or  reaction  rate  flux.  However,  the  odd  parity 
equation  and  flux  has  the  odd  information  which  is  impor¬ 
tant  in  determining  the  net  current. 

In  order  to  see  that  this  is  indeed  the  case  the 
reader  can  easily  verify  that  in  any  orthogonal  coordinate 
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system  the  surface  term  is  an  odd  function  on  the  unit 
sphere.  It  is  an  odd  function  in  the  direction  or 
equivalently  in  the  angles  p  and  x*  Then,  by  writing  the 
angular  flux  as  in  Eq  (84)  and  the  parity  flux  as  in  Eqs 
(29)  and  (30),  he  can  show  (Wills,  1981)  that 

/Vu( r,Q)dr  =  0  (88) 

•'VTT 

where 

A(r,fi)d6  =  A(r,-fbd6  (89) 

Jn it 

and  the  integration  is  carried  out  over  all  directions 
(4-rr  steradians)  .  Therefore 

=  /,'Fg(  r,6)d6  (90) 

•/Mr  .A  TT 

Furthermore,  since  the  term  n-n  is  represented  by 
odd  spherical  harmonic  functions  (Wills,  1981)  we  can  write 


f  (6.fi)¥8(r,n)d6  =  0 

J  4TT 


(91) 
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and  the  particle  current  as 


f  (6-n)4>(r  f«)dfl  =  f  (Q*n)Yu(r  t8)dn  (92) 

J'n r 


A  A  ^  A  A 

where  we  used  the  fact  that  Q*n  and  ¥  (r,^)  are  odd,  whereas 
^(r,fi)  is  even  on  the  unit  sphere. 

This  odd-even  property  of  the  second  order  formulation 
is  also  important  in  the  choice  of  a  finite  element  trial 
basis.  We  propose  to  use  a  dual  basis  --  a  piecewise  poly¬ 
nomial  (spline)  basis  in  the  spatial  variables  and  spherical 
harmonics  in  angle.  Our  choice  of  spherical  harmonics,  or 
an  equivalent  tensor  product  of  surface  harmonics  (Sansone, 
1977)  in  the  angular  variables  is  based  on  the  following: 

i.  The  odd-even  properties  of  the  second  order  flux 
and  equations . 

ii.  The  usual  practice  of  expanding  the  scattering 
cross-sections  in  spherical  harmonics  (Bell 
and  Glasstone,  1970). 

iii.  The  streaming  and  surface  normal  terms  are 
expanded  by  the  use  of  Legendre  functions. 


iv.  The  orthogonal  and  parity  characteristic 
of  spherical  and  surface  harmonics, 
v.  Spherical  harmonic  functions  have  been 

successful  in  mitigating  ray  effects  (Reed, 
1972). 

We  also  chose  a  spatial  finite  element  basis  in  order 
to  ease  the  computational  burden  of  the  five-dimensional 
problem.  A  finite  element  basis  can  easily  model  diffi¬ 
cult  and  complicated  geometries.  Furthermore,  a  finite 
element  spline  basis  is  easy  to  construct  and  their 
piecewise  nature  facilitates  the  evaluation  of  the  many 
integrals  in  the  Galerkin  weak  form.  There  is  also  a 
wide  selection  of  spline  basis  functions  and  element  shapes 
to  choose  from  (Prenter,  1975).  This  provides  the 
flexibility  which  is  required  in  solving  general  particle 
transport  problems . 

Expanding  our  solution  in  this  dual  basis  of  poly¬ 
nomial  splines  and  spherical  harmonics  we  can  now  write 
the  trial  function  expansion  of  Eq  (65)  and  Eq  (87)  with 


Nt(f,fi)  =  (93) 

where  -Usms+Jl  and  the  index  i  varies  with  Z ,  m  and  j. 
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This  is  a  global  basis  of  harmonic  functions  in  angle  and 
a  local  finite  element  basis  in  the  spatial  variables  which 
represent  r.  Referring  to  Figure  3,  these  could  be  the 
x,  y,  z  variables  in  rectangular  geometry  where  the  angle 
or  direction  variables  are  u  and  x*  Furthermore,  if  there 
is  symmetry  in  the  problem  we  could  eliminate  some  of  the 
independent  variables  and  thus  reduce  the  dimension  of  the 
trial  space. 

Spherical  Harmonics  and  Surface  Harmonic  Tensors. 
Spherical  harmonics  are  complex  functions  in  the  direction 
variables  of  fi.  These  variables  are  usually  denoted  by 
the  symbols  p  and  x-  In  rectangular  geometry  they  are  de¬ 
fined  as  in  Figure  3.  Their  definition  in  other  orthogonal 
geometries  can  be  found  in  Bell  and  Glasstone  (1970).  They 
specify  directions  or  points  on  a  unit  sphere  where  ue[-l,  1] 
and  xe[0,  2tt],  The  normalized  spherical  harmonic  functions 
are  defined  as  (Wills,  1981) 


+  iQ 


s 

l  ,m 


(94) 
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where  P^m^u)  are  the  associate  Legendre  functions 
(MacRobert,  1967)  and 


I 

r  -  lil+l  q-m)  I 

c)ifm  y-znr  *  q+m) : 


(95) 


We  have  also  defined  Q^m  and  Q^m  to  be  the  even  and  odd 
surface  harmonic  tensors 


<4,  *  ctaPHm(u)cos(mX) 


(96) 


and 


C£mPim(u)8ln<nX) 


(97) 


The  usual  spherical  harmonic  expansion  of  a  function 


f(ft)  is  written  as 


(98) 


f<6) 


L  +Z 

-  I  I  WWfi> 

£=0  nr=-£ 


where  f  ^  m  are  the  expansion  coefficients  which  are 
determined  by  the  inner  product 


£»  *  <£<a> ^  jff(Sn*w(S)dsi 


and  Y*m  is  the  complex  conjugate.  We  can  rewrite  Eq 
without  a  loss  of  generality  (Wills,  1981)  as 


f(fi) 


i  i  nA  *  ££°mc 

i= 0  m=0 


A  O 

where  f^m  +  f^m  and  we  have  used  the  fact  that 


l,  -m 


(-l)mY 


Jim 


=  ( -l)mQC 


Urn 


+  K-l)m+1Q?m 


tjfa 


(99) 
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(100) 
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Using  a  spherical  harmonic  expansion  in  Eqs  (65) 
and  (87)  produces  a  matrix  problem  which  is  complex  if 
there  is  no  symmetry  in  the  angle  x  (Blomquist  and  Lewis , 
1980).  The  matrix  of  the  Galerkin  weak  form  is  then  a 
Hermetian  matrix  and  the  computations  must  be  carried  out 
in  complex  arithmetic.  However,  since  we  know  that  the 
sources  are  all  real  and  that  the  solution  is  real,  we 
could  use  an  alternate  and  equivalent  real  trial  space 
(Fletcher,  1983).  This  trial  space  is  just  the  surface 
harmonic  tensors  of  Eq  (94)  where  the  expansion  of  Eq 
(100)  is  now  written  as 


f  (Q) 


"S  + 


s  s 

£mxi!m  ' 


(102) 


and  the  trial  function  of  Eq  (93)  as 


Ni(r,fi)  =  Q^m  +  Bj(r)  with  0<m<«-  (102a) 


In  order  to  see  that  the  expansions  of  Eqs  (98)  and 
(102)  are  equivalent,  we  need  to  separate  the  spherical 
harmonic  coefficients  into  their  real  and  imaginary  parts 
(Courant  and  Hilbert,  1953)  as 
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for  m>0 
for  m<0 


(103) 


By  substituting  Eq  (103)  into  Eq  (98)  we  can  derive  the 
real  expansion,  Eq  (102).  These  harmonic  functions  are 
complete  and  converge  uniformly  and  in  the  mean  (Churchill 
and  Brown,  1978).  They  are  also  orthogonal  in  that  for 
A  f  A1  and  m  #  m' 


<YAm<6>*Vm^»  - 


(B  <0C  0S  > 


"  <(4X’m’> 


=  0 


(104) 


Another  important  property  of  these  harmonic  func¬ 
tions  is  that  they  have  special  parity  characteristics. 

This  is  a  very  important  property  computationally  and  in 
choosing  the  angular  trial  space.  It  allows  us  to  construct 
a  basis  which  is  best  suited  to  accomodate  the  odd- 
even  (parity)  properties  of  the  even  and  odd  parity 
fluxes.  In  solving  the  even  or  odd  parity  equations  we 
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need  to  use  a  trial  space  which  is  complete  in  either  the 
even  or  odd  functions  on  the  unit  sphere.  In  effect,  we 
need  to  limit  our  trial  space  to  either  even  or  odd  func¬ 
tions.  This  restriction  is  implicit  in  our  definition  of 
the  parity  flux  (Eqs  (29)  and  (30)).  There  we  separated 
the  angular  flux  into  a  sum  of  components  (the  even  and 
odd  parity  fluxes)  which  are  even  and  odd  on  the  unit 
sphere . 

Therefore,  in  solving  Eq  (87),  we  need  to  use  an  an¬ 
gular  trial  space  of  even  functions.  By  our  definition, 
the  odd  components  in  this  function  space  are  all  zero. 
Using  the  parity  characteristic  of  spherical  harmonics 
this  trial  space  is  just  the  expansion  of  Eqs  (98)  or  (102) 
where  the  index  £  is  restricted  to  be  even  (Davis,  1966). 
For  a  solution  to  the  odd  parity  problem  we  could  construct 
a  similar  reduced  trial  space  where  £  is  odd.  In  either 
case  we  have  reduced  our  trial  space  to  a  subset  of  the 
larger  solution  space  which  is  required  in  the  Pn  method 
(Greenspan  et  al.,  1968).  Furthermore,  this  smaller  trial 
space  reduces  the  size  of  the  matrix  problem  which  must  be 
solved.  This  in  turn  reduces  the  computational  effort 
which  is  required. 

In  order  to  demonstrate  that  these  harmonic  functions 
have  parity  we  need  to  show  that 


(105) 


or  equivalently  that 


C<-S>  - 


(106) 


and 


Qsta(-S)  -  <-l)eQ=m<S)  <107> 


Here,  specifies  the  direction  opposite  to  6  and  can  be 
represented  by  -y  and  it  +  x  (See  Figure  3).  The  proof  of 
these  relationships  is  straightforward.  It  includes  the 
fact  (Hobson,  1965)  that 


(108) 


102 


(109) 


eim(*  +  X)  =  (_i)meimx 


sin[m(7T  +  x)l  =  (-l)msin(mx)  (110a) 


and 


cos[m(Tr  +  x)J  =  (-l)mcos  (mx)  (110b) 


Finally,  we  can  show  that  in  the  limit  of  diffusion 
theory,  the  Galerkin  second  order  equations  reduce  to  the 
diffusion  equation  and  Fick's  law.  This  equivalence  is 
important  in  that  it  shows  the  relationship  between  our 
finite  element  projection  method  and  the  P-^  method.  It 
also  shows  that  our  approach,  in  the  limit  of  a  linearly 
anisotropic  approximation,  is  consistent  with  these 
established  and  proven  techniques .  This  equivalence  is 
presented  and  discussed  in  Appendix  D.  Furthermore,  we  can 
also  expand  the  streaming  operators  and  cross-sections  of 
the  second  order  forms  by  using  the  well-known  recursive 
properties  of  spherical  harmonics.  These  relationships 
are  included  in  Appendix  E. 


The  Muitigroup  Solution.  A  complete  solution  of 
general  particle  transport  problems  must  include  a 
rigorous  treatment  of  the  time  and  energy  dependences . 

As  shown  in  Chapter  IV,  the  usual  approach  is  simply  to 
extend  the  steady-state  and  time -independent  solution  to 
energy  and  time  dependent  problems .  For  the  second  order 
transport  equations  this  includes  the  usual  multigroup 
approach  in  energy  and  a  weighted  difference  or  Crank - 
Nicholson  scheme  in  time  (Duderstadt  and  Martin,  1979). 

In  order  to  keep  the  algebra  and  notation  to  a 
minimum  and  at  a  manageable  level  we  will  not  rewrite 
the  Galerkin  weak  form  with  an  explicit  time  and  energy 
dependent  notation  (subscript).  Instead,  we  will  refer 
the  reader  to  Duderstadt  and  Martin  (1979)  and  to 
Chapter  IV.  However,  we  will  discuss  a  time  and  energy 
dependent  solution  within  the  context  of  Eq  (87).  Our 
intent  is  to  discuss  the  important  aspects  of  a  general 
solution  which  includes  the  usual  time  and  energy  depen¬ 
dent  treatment. 

The  discrete  time  and  energy  dependent  equations  are 
just  the  system  of  Eq  (87)  with  redefined  operators,  cross- 
sections  and  source  terms  (See  Chapter  IV) .  The  even  parity 
flux  is  now  defined  with  an  explicit  time  and  energy  (multi- 


group)  dependence.  The  source  terms  will  include  known 
flux  values  from  a  previous  time  step  and  sources  from 
other  energy  groups.  This  problem  will  be  implicit  in  time 
with  the  solutions  at  a  given  time  step  being  dependent  on 
those  at  previous  times.  Furthermore,  the  multigroup 
equations  are  coupled  by  the  source  terms  whereby  the 
sources  in  a  given  group  are  dependent  on  the  flux  in 
other  groups. 

The  general  solution  strategy  is  to  begin  with  some 
initial  condition  and  to  proceed  by  a  marching  scheme  in 
time.  At  each  time  step  we  then  solve  the  multigroup 
energy  dependent  problem  by  the  usual  iterative  procedure. 
This  multigroup  problem  is  a  set  of  coupled  one  group 
problems  where  the  parity  sources  of  Eq  (87)  are  energy 
and  time  dependent  sources.  These  sources  can  be  written 
in  a  form  similar  to  Eq  (26)  as 


Qg  =  sg  +  sg  +  sg  +  KUt,yg) 
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Qv  -  Sj  +  SU  +  sg  +  K(At,^) 


[i. 


Here  we  have  omitted  the  explicit  time  index  and  included 
fissions  and  group  scatter  sources.  KCVt,^)  and  K(Vt,’l'^) 
represent  terms  which  depend  on  the  time  discretization 
that  is  used  and  b  is  the  energy  group  index.  ^  and 
are  the  even  and  odd  parity  fluxes  from  the  previous  time 
step  or  initial  conditions,  and 

»  even  parity  downscatter  source. 

S®  s  even  parity  upscatter  source, 
g 

■  even  parity  fission  source. 
s|j  =  odd  parity  downscatter  source. 
sJJ  **  odd  parity  upscatter  source. 

*  even  parity  inhomogenous  source. 

*  odd  parity  inhomogenous  source. 

These  sources  are  similar  to  the  usual  multigroup 
sources  where,  for  example,  the  odd  parity  downscatter 
source  for  energy  group  b  can  be  written  as 


where  l^biG  and  group  one  has  the  most  energetic  particles. 
The  other  sources  can  be  similarly  defined  where  the 
scattering  cross-sections  are  defined  by  Eqs  (33)  and  (34). 

It  is  important  to  note  that  because  fissions  occur 
isotropically  they  only  contribute  to  the  even  parity 
source  term.  Furthermore,  with  isotropic  inhomogenous 
sources  and  scattering  the  odd  parity  sources  are  all  zero. 
However,  for  problems  with  anisotropic  scattering  we  must 
include  and  compute  the  odd  parity  downscatter  and  up- 
scatter  sources.  These  and  the  K(Vt,^)  term  of  Eq  (112) 
are  expressed  in  terms  of  the  odd  parity  flux.  Therefore  we 
must  compute  the  odd  parity  group  fluxes  in  order  to 
determine  the  source  terms  for  anisotropic  and  time  depen¬ 
dent  problems. 

This  is  a  disadvantage  of  our  second  order  formula¬ 
tion.  If  we  are  only  interested  in  scalar  flux  values  and 
thus  a  solution  to  the  even  parity  equation  then  for 
anisotropic  and  time  dependent  problems,  we  must  indirectly 
compute  the  odd  parity  flux.  This  can  be  accomplished  by 
inserting  Eq  (37)  into  Eqs  (112)  and  (113).  However,  this 
increases  the  computational  work,  especially  for  problems 
with  upscatter  sources.  Nonetheless,  we  are  not  computing 
4^  explicitly,  and  a  solution  can  be  obtained  without  having 
to  solve  the  odd  parity  equations. 


Except  for  the  extra  effort  which  is  required  in 
computing  the  multigroup  and  time  dependent  sources  the 
solution  approach  is  similar  to  that  which  is  normally  used 
to  solve  the  first  order  transport  equation  (Sanchez  and 
McCormick,  1982).  This  involves  an  inner  and  outer  itera¬ 
tion  strategy.  At  each  time  step  we  carry  out  a  multigroup 
solution  where  we  solve  the  one  group  steady  state  matrix 
problem  for  each  energy  group.  This  is  the  inner  iteration 
where  we  calculate  the  space  and  angle  flux  distribution 
T^(r,fi)  for  each  energy  group  b.  The  usual  practice  for 
solving  the  diffusion  equation  is  to  use  an  iterative 
(indirect)  solver  to  solve  the  matrix  problem  (Bell  and 
Glasstone,  1970).  This  practice,  together  with  the  itera¬ 
tive  procedure  of  the  discrete  ordinate  method,  has  led  to 
the  concept  of  an  "inner  iteration  strategy." 

An  outer  iteration  provides  a  solution  for  the  group 
fluxes  ^(r,^)  for  all  of  the  energy  groups.  It  involves 
a  sweep  through  all  energy  groups  in  order  of  decreasing 
energy  to  provide  an  updated  steady  state  solution.  The 
outer  iteration  is  then  repeated  until  some  convergence 
criteria  is  satisfied  or  a  prescribed  number  of  iterations 
are  performed.  For  problems  which  do  not  have  fissions 
or  upscatter  sources,  only  one  outer  iteration  is  required. 
The  group  flux  T®(r,fi)  converges  in  one  outer  iteration. 
However,  fissions  and  upscatter  sources  depend  on  the  un- 
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known  and  to  be  computed  flux  in  lower  energy  groups. 
Therefore,  the  usual  practice  is  to  begin  with  an  initial 
flux  guess  for  all  of  the  energy  groups  and  to  perform  a 
number  of  outer  iterations  before  convergence  is  achieved. 
At  each  outer  iteration  the  fissions  and  upscatter  sources 
are  computed  by  using  the  flux  solutions  of  the  previous 
outer  iteration. 

The  Matrix  Problem.  A  numerical  solution  of  the 
second  order  transport  problem  requires  that  we  solve  a 
matrix  problem.  This  matrix  problem  is  obtained  by  expan¬ 
ding,  evaluating  and  assembling  the  individual  terms  of 
the  Galerkin  weak  form,  Eq  (87).  This  involves  the 
evaluation  of  a  number  of  integrals  and  the  construction 
of  a  problem  matrix  and  source  vector.  Finally  we  must 
solve  this  system  for  the  flux  expansion  coefficients 
(<f>^).  These  coefficients  can  then  be  substituted  into  Eq 
(65)  to  give  the  even  parity  flux. 

The  matrix  problem  can  be  written  as 


[a]  [♦]  -  [s 


(114) 


where 


A  =  positive  definite  symmetric  problem 
matrix 


$  =  flux  expansion  coefficient  vector 
S  =  right  hand  side  or  source  vector 

To  get  from  Eq  (87)  to  (114)  we  will  follow  the  usual 
finite  element  practice  (Huebner,  1975)  of  evaluating  the 
integrals  by  Gauss -Legendre  quadrature.  Numerical  inte¬ 
gration  with  Gauss  quadrature  can  be  written  as  (Atkinson, 
1978) 


fa f(x)dx=£  w.f(x. 

i=l  J  J 


(115) 


where  the  w^ ' s  are  the  Gauss  Legendre  weights  and  x^ ' s 
are  the  zeros  of  the  nth  degree  Legendre  polynomial.  The 
quadrature  formulae  Eq  (115)  integrates  all  polynomials  of 
degree  <2n-l  exactly. 

In  order  to  reduce  the  computational  effort  a  complete 
expansion  of  the  Galerkin  weak  form  and  identification  of 
the  distinct  integrals  is  required  (Wills,  1981).  These 
integrals  are  evaluated  once  and  then  used  to  assemble  the 
problem  matrix  and  source  vector.  We  can  then  use  an 
indirect  solver  (Hageman  and  Young,  1981),  or  Cholesky 
decomposition  (Stewart,  1973)  to  solve  this  positive  de- 
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finite  symmetric  system.  Furthermore,  the  problem  matrix 
will  have  a  large  number  of  zeros  due  to  the  use  of  finite 
elements.  Therefore,  we  could  use  a  sparse  matrix  storage 
scheme  (Jennings,  1977)  to  reduce  the  problem  size.  This, 
coupled  with  the  fact  that  the  problem  matrix  is  symmetric, 
will  decrease  the  computer  storage  requirements  and  costs. 


CHAPTER  VIII 


THE  AIR-OVER-GROUND  PROBLEM 


We  now  apply  our  method  to  a  real  and  difficult  pro¬ 
blem  —  the  air-over-ground  problem.  This  problem  involves 
the  deep  penetration  of  particles  in  the  atmosphere  along 
with  reflection  and  propagation  into  the  ground.  These  par¬ 
ticles  are  neutrons,  gamma  rays  and  x-rays.  They  propagate 
away  from  a  point  source  in  a  two-dimensional  air-over¬ 
ground  geometry.  Solutions  to  this  problem  have  a  variety 
of  applications.  These  applications  include  the  prediction 
of  collateral  damage  and  radiation  exposure  in  a  nuclear 
environment . 

The  air-over-ground  problem  increases  in  complexity 
because  it  has  point  sources,  forward  peaked  anisotropic 
scattering,  an  exponentially  varying  atmosphere,  and  an  air- 
ground  interface  (Pace  et  al.,  1975).  Numerical  solutions 
to  this  problem  already  exist.  The  main  solution  techniques 
are  Monte  Carlo  and  discrete  ordinates.  However,  both 
Monte  Carlo  and  discrete  ordinates  have  severe  difficulties 
and  limitations  (Loewe  et  al.,  1983).  A  Monte  Carlo  solu¬ 
tion  has  statistical  errors.  These  are  expected  uncertain- 


ties  in  the  approximate  solution.  However,  these  uncertain¬ 
ties  become  large  as  the  problem  dimensions  increase. 
Therefore,  over  large  space  regions  the  Monte  Carlo  solution 
has  poor  statistics  and  unacceptable  errors.  It  also  re¬ 
quires  many  hours  of  computer  execution  time.  Discrete 
ordinates  is  also  limited  by  deficiencies  in  the  difference 
scheme,  ray  effects  and  long  execution  times. 

Because  of  these  difficulties  many  researchers  have 
investigated  and  developed  alternate  solution  techniques 
for  solving  the  air-over-ground  problem.  Roberds  and 
Bridgman  (1977)  have  used  a  projection  of  "specially 
tailored"  trial  functions  to  solve  this  problem.  Sbulstad 
(1976),  Eamon  (1976)  and  Sacenti  and  Jacobsen  (1975)  have 
used  a  mass  integral  scaling  technique.  In  this  approach 
scale  factors  are  used  to  correct  infinite  homogenous 
air  results.  Souders  (1981)  has  also  applied  projection 
methods  to  the  even  parity  transport  equation.  He  attempted 
a  collocation  and  Galerkin  solution  and  concluded  that 
these  methods  could  not  be  successfully  used  to  solve  this 
problem.  In  contrast  to  our  use  of  spherical  harmonics  he 
used  Lagrange  polynomials  and  special  ellipsoid  trial 
functions  to  model  the  angular  dependences.  He  did  not 
include  the  air-ground  interface  or  a  nonhomogenous 
atmosphere. 
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The  most  general  solution  of  the  air-over-ground  pro¬ 
blem  can  be  obtained  in  a  three  dimensional  spatial  geome¬ 
try.  This  is  a  time  and  energy  dependent  problem  with  two 
angle  variables  (u,x).  Therefore,  the  general  air-over¬ 
ground  problem  requires  a  detailed  seven  dimensional 
solution.  However,  if  in  cylindrical  geometry,  we  confine 
the  point  sources  to  the  z-axis ,  then  the  problem  becomes 
two  dimensional  in  the  spatial  variables.  It  has  azimu¬ 
thal  symmetry  and  can  be  solved  as  a  six  dimensional  problem 
in  a  cylindrical  (r,z)  geometry. 

It  is  not  our  purpose  to  provide  the  most  general 
solution  to  the  air-over-ground  problem  here.  Instead,  we 
will  reduce  the  problem  further  to  a  monoenergetic  steady- 
state  problem  in  cylindrical  geometry.  The  point  sources 
will  be  confined  to  the  z-axis.  However,  we  will  include 
first  scatter  sources,  anisotropic  scattering,  an  exponen¬ 
tially  varying  atmosphere  and  the  air-ground  interface. 
Furthermore,  our  solution  can  be  easily  extended  to  the 
general  time  and  energy  dependent  (multigroup)  problem. 

As  was  pointed  out  in  Chapter  VII,  this  is  a  straight¬ 
forward  extension  which  requires  the  solution  of  a  number 
of  coupled  steady-state,  monoenergetic  problems.  Therefore, 
our  primary  objective  is  to  demonstrate  that  the  finite 
element  projection  method  can  be  used  to  solve  the  steady- 
state  one  group  air-over-ground  problem. 
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Problem  Geometry  and  Boundary  Conditions 


This  problem  will  be  modeled  in  cylindrical  (r,z) 
geometry  with  vacuum  boundary  conditions.  It  is  a  steady- 
state,  monoenergetic  (one  group)  four -dimensional  problem 
with  two  spatial  (r,z)  and  two  angle  (u,x)  variables.  The 
coordinate  system  and  directional  angles  of  this  geometry 
and  the  air-over-ground  problem  are  shown  in  Figure  (4). 
y  is  the  cosine  of  the  angle  formed  by  the  z-axis  and  the 
particle  velocity  vector  .  X  is  the  angle  between  the 
planes  formed  by  the  r  vector  and  z-axis  and  that  of  the 
vector  and  z-axis. 

For  the  air -over -ground  problem  the  scattering  cross- 
sections  will  be  highly  peaked  in  the  forward  direction 
(Bartine  et  al. ,  1977).  This  is  due  to  the  scattering  pro¬ 
perties  of  air  and  the  high  particle  energies  which  exist 
in  this  problem.  Because  of  this  the  exterior  boundary 
condition  for  this  problem  will  be  approximated  by  a 
vacuum  boundary  condition  (Burgio,  1975).  This  boundary 
condition  is  represented  by  Eq  (9)  and  for  the  second  order 
transport  equation  by  Eqs  (46)  and  (47).  In  the  Galerkin 
weak  form  and  the  equivalent  variational  problem  it  is 
defined  by  Eq  (64a). 

Furthermore,  in  two-dimensional  cylindrical  (r,z) 


geometry  there  is  symmetry  in  the  angle  x-  This  symmetry 
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can  be  written  as 


<J>(r,fi  )  =  4>(r,6  )  for  Q 

1  2  p  2 


Pi 


2(n-fi  )n 
Pi 


(116 


or  as 


<}>  (r  ,u  ,x)  -  <J>(r,y,-x)  -  $(r,y,2ir-x)  (117 


The  symmetry  x  is  shown  in  Figure  5 ,  where  the  vectors  n 
and  p  are  perpendicular. 

Note  that  because  of  the  problem  geometry  and  the 
exponentially  varying  air  density  (in  the  z  direction) , 
only  azimuthal  symmetry  in  the  angle  x  is  assured.  There 

A  A 

is  no  symmetry  in  y(cos9)  and  therefore  <J>(r,f2)  will  not 
be  equal  to  <J>(r,-J2).  The  symmetry  condition,  Eqs  (116) 
and  (117),  implies  that 


<J>(r,y,x)  =  even  functions  in  x 
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A  Schematic 


of  the  Angular  Symmetry. 
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The  problem  geometry  and  coordinate  system  as  shown 
in  Figure  A ,  together  with  the  symmetry  in  x  requires  that 
when  P  is  equal  to  zero  the  solution  is  a  constant  in  X- 
This  also  implies  that  at  p  =  0  the  angle  x  must  be  equal 
to  zero.  Therefore,  along  the  z-axis  we  have  that 

4>(r,Q)  =  0  for  o  =  0  (120) 

and  x  =  0 

At  the  air-ground  interface  the  angular  flux  is 
continuous  in  the  spatial  variable  r.  This  interface 
condition  which  is  represented  by  Eq  (13)  is  a  natural 
condition  of  the  Galerkin  weak  form  (See  Appendix  C) . 
However,  this  is  a  plane  interface  for  which  it  can  be 
shown  (Bell  and  Glasstone,  1970)  that  the  angular  flux  is 
discontinuous  in  the  variable  y.  This  discontinuity  exists 
at  y  =  0  and  depends  primarily  upon  the  problem  geometry, 
sources  and  material  properties.  The  most  severe  cases 
are  those  of  a  vacuum  boundary  and  interfaces  where  the 
adjoining  materials  are  very  dissimilar  in  terms  of 
sources  and  cross-section  data.  The  air-over-ground  pro¬ 
blem  is  in  this  later  category. 


It  is  important  to  note  that  this  discontinuity  does 
not  present  a  difficulty  in  the  Galerkin  weak  form.  This 
formulation  treats  the  interface  condition  as  a  natural 
condition  and  would  provide  an  approximate  solution  to  the 
problem.  However,  our  trial  space  must  be  complete  and 
in  this  case  completeness  includes  the  ability  to  approxi¬ 
mate  discontinuous  functions.  For  problems  where  this 
discontinuity  is  severe  the  choice  of  a  reasonable  and 
"admissible"  trial  space  is  important. 

Our  trial  space  of  harmonic  functions  is  complete. 

The  associated  Legendre  polynomials  are  also  complete. 

They  can  approximate  discontinuous  functions  in  the  mean  by 
an  infinite  series  where  the  number  of  functions  (N)  tend  to 
°°  (Bell  and  Glasstone,  1970).  However,  we  can  get  a  better 
approximation  with  a  smaller  trial  space  (N) ,  if  we  include 
functions  which  are  discontinuous  at  y  =  0.  This  is  the 
basic  idea  of  the  Double  Pn  approximation  (Duderstadt  and 
Martin,  1979).  We  simply  separate  the  trial  space  in  two 
half  ranges  in  y  and  use  separate  Legendre  expansions  on  each 
range,  i.e.,  ye [-1,0)  and  ye [0,1].  With  this  approach  we 
can  easily  treat  problems  with  severe  interface  discon¬ 
tinuities.  A  very  accurate  solution  of  the  air-over-ground 
problem  may  require  such  a  treatment. 


The  Finite  Element  Basis  and  Problem  Equations 

Applying  the  finite  element  projection  method  of 
Chapter  VII  to  the  air -over -ground  problem  requires  a  nu¬ 
merical  (computer)  solution  of  the  Galerkin  weak  form  Eq 
(87).  However,  we  must  further  define  and  choose  our 
trial  space  in  accordance  with  the  problem  symmetries 
and  boundary  conditions.  In  this  solution  we  will  follow 
the  usual  practice  of  not  using  a  Double  Pn  approximation 
in  angle.  This  approximation  may  only  be  needed  for  those 
problems  where  a  very  accurate  solution  at  the  air-ground 
interface  is  desired. 

The  general  solution  strategy  of  Chapter  VII  employs 
a  dual  basis  of  global  harmonic  functions  in  angle  and 
finite  elements  in  the  spatial  variables.  We  could  also 
use  four  dimensional  finite  elements  in  both  space  and  an¬ 
gle  (Miller  et  al. ,  1973).  However,  these  four  dimensional 
phase-space  finite  elements  will  be  very  costly  (computer 
costs)  and  inefficient  (Wheaton,  1978).  This  is  due  to 
the  added  complexity  of  anisotropic  scattering.  Aniso¬ 
tropic  scattering  increases  the  data  management  and 
computational  difficulties.  A  local  basis  in  angle  requires 
that  the  scattering  contribution  to  each  element  must  be 
computed  on  an  element  by  element  basis  for  all  space  and 
angle  elements  within  the  problem  domain.  Therefore,  a 
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four-dimensional  phase-space  finite  element  formulation  of 
this  problem  is  not  a  very  attractive  or  realistic  approach. 

We  will  use  a  dual  basis  of  spherical  harmonics  and 
splines.  The  global  nature  of  the  spherical  harmonic 
functions  would  facilitate  the  treatment  of  anisotropic 
scattering.  This  expansion  which  must  also  satisfy  the 
symmetry  conditions  of  Eqs  (117)  and  (120)  can  be  written 
as 


^(z.p.y.x)  - 


IR  L  A 

ISZ  A1Blz(i!)Blr(p)Q 
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where  A.  are  the  expansion  coefficients  and  B.  (z)B.  (p) 
form  a  tensor  product  of  normalized  B-splines  on  a  rectan¬ 
gular  (P,z)  grid  (Deboor,  1978).  Q^m  is  a  surface  harmonic 
function  which  is  defined  by  Eq  (96)  with  SL  restricted  to 
be  even. 

The  definition  of  a  polynomial  B-spline  is 


B(x)  = 
i,k,t 


(ti4k"  ti)[ti' 


ti+klGk(x’t) 


(122) 


where 


Gk(x,t) 


(t 


k-1 


(t-x)^  ■*"  for  t  i  x 
0  for  t  £  x 
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and 

k  =  1,2....=  order  of  the  B-spline 
i  =  index  of  the  i-th  normalized  B-spline 

t^  and  t^+k  represent  the  node  points  and  [t^ . fci+k^  <*e“ 

fines  the  usual  k-th  Newton  divided  difference.  The  graph 
of  an  unnormalized  cubic  spline  function  with  evenly 
spaced  nodes  (knots)  is  shown  in  Figure  6.  High  order 
splines  (k  z  2)  are  (k-2)  times  continuously  differentiable. 
These  are  very  smooth  functions,  however,  we  can  also 
construct  low  order  linear  splines  and  splines  (k  =  1) 
with  jump  discontinuities  at  the  nodes  (Schultz,  1973). 

B-splines  have  been  used  extensively  as  interpolation 
functions  (Schumaker,  1981).  Spline  interpolation  requires 
the  solution  of  a  smaller  system  of  equations  than 
interpolation  with  Hermite  or  Lagrange  polynomials.  The 
usual  practice  is  to  use  normalized  B-splines  in  a  N- 


dimensional  space  whereby 


N 

y,  B.(g)^  =  1  for  all  x 


(124) 


A  tensor  product  of  B-splines  is  being  used  here  for 
the  following  reasons : 

1.  They  are  piecewise  continuous  and  form  a  local 
basis  such  that  the  integral  B^(x)B^ (x)  dx  is  zero 
if  |i  -  j|>k  where  k  is  the  order  of  the  B-splines. 

This  reduces  the  number  of  integrals  which  must 

be  evaluated  and  also  produces  a  sparse  and  banded 
coefficient  matrix. 

2 .  A  separation  of  the  p  and  z  integration  variables 
is  possible. 

3.  For  a  given  problem  partition  (mesh  spacing)  poly¬ 
nomial  splines  will  produce  a  coefficient  (problem) 
matrix  that  is  smaller  but  less  sparse  than  Hermites 
or  Lagrange  polynomials . 

In  the  Galerkin  weak  form  of  Eq  (87)  the  test  or  weight 
functions  can  now  be  defined  as 
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N^z.p.vux)  *  Bjz(z)B.r(p)(^n 


(125) 


The  spatial  trial  functions  of  B-splines  are  defined  in 
the  global  problem  coordinates .  In  this  solution  a  local 
element  coordinate  system  and  a  parametric  reoresentation 
will  not  be  used.  Instead  we  will  use  ring  type  axisymme- 
tric  rectangular  elements  which  are  defined  in  the  problem 
coordinates.  These  are  four  node  two-dimensional  elements 
with  axial  symmetry.  They  will  be  constructed  on  a  non- 
uniform  mesh  in  the  (P,z)  plane.  Furthermore,  in  this 
representation  the  expansion  coefficients  of  Eq  (121)  do  not 
represent  the  solution  values  at  the  nodes  except  when 
linear  B-splines  (k  =  2)  are  used.  In  the  terminology  of 
Chapter  III  and  with  this  exception,  these  coefficients 
are  generalized  coordinates . 

By  substituting  Eqs  (125),  (121)  and  (64a)  into  Eq 
(87)  we  get 


IZ  IR 

y  i 

&lir-l 


i[/r{<S-7(V  <&> 


+  ^kn^'VO^  +  W 


n’n'BjQknBiQ£mdfids 


-/R{<5-7<Vkn  >'kU«U>  +  <Bj(&’<38>}d* 


(126) 


Here  Z  and  k  are  even  integers  and  for  simplicity  the  p,  z, 
y,  X  dependences  are  omitted  and 


B.  = 

i 


Bi2<z)Bir(P> 


(127) 


B 


j 


(z)Bjr(p) 


(128) 


Eq  (126)  can  be  written  as  a  matrix  problem  in  the 
form  of  Eq  (114)  where  the  A^s  are  the  flux  expansion 
coefficients.  The  elements  of  the  problem  matrix  and 
source  vector  are  obtained  by  evaluating  and  summing  the 
individual  expanded  terms  of  Eq  (126).  Details  of  this  ex¬ 
pansion  are  carried  out  in  Appendix  F  where  the  surface 
normal  and  (streaming)  gradient  operators  are  defined  in 
cylindrical  geometries. 

The  directional  gradient  operator  in  cylindrical 
geometry  is  defined  as  (Bell  and  Glasstone,  1970) 


n*v<j> 


cos (x) 


9  (p<p )  _  1^  9<  4»  Vl-y 2sinx  f  +  y  9<J> 
9  p  "  p  9X  Tz 


(129) 
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This  is  the  conservative  form  of  the  directional  derivative 
in  two  dimensional  (p, z)  geometry  with  azimuthal  symmetry. 
The  Ku  and  G®  operators  have  been  defined  in  Eqs  (42)  and 
(39).  The  scalar  product  of  the  velocity  and  normal  vec¬ 
tors  can  be  written  as 

|fl*n  for  the  horizontal  outer  surfaces  (top  or 
z  bottom)  of  the  problem  cylinder,  and  as 

ft*n  for  the  vertical  surface  (side)  of  the 

p  cylinder  (130 

where 

y  on  the  top  surface,  and 

A  A 

n*n  “ 

-y  on  the  bottom  surface  (131 

and 

fi-np  =  Vl-y2  cosx  (132 

The  normal  unit  vectors  ri  and  h  are  shown  in  Figure  7 . 

p  z  ° 

Expanding  the  expressions  in  Eq  (126)  produces  an  in¬ 
tegral-differential  equation  which  has  twenty-eight  terms 
(See  Appendix  F) .  These  terms,  except  for  the  source  terms 

can  easily  be  separated  into  a  product  of  z,  p,  y  and  x 


integrals.  This  is  an  integral  separation  of  variables 
which  is  a  direct  result  of  Eq  (121);  where  it  is  assumed 
that  the  solution  can  be  expressed  in  a  form  where  the 
dependent  variables  are  separable.  This  separation  pro¬ 
perty  simplifies  the  individual  integrals  which  have  to 
be  evaluated.  It  allows  for  the  evaluation  of  only  single 
integrals  and  not  the  more  complicated  double,  triple  or 
quadruple  integrals.  By  this  separation  of  variables  it 
may  be  possible  to  integrate  most  of  these  single  integrals 
analytically  and  thus  avoid  a  numerical  integration  process. 

The  term  by  term  expansion  of  Eq  (126)  has  produced 
thirty-seven  distinct  single  integrals.  These  integrals 
can  be  found  in  Appendix  G.  The  angle  integrals  are  only 
dependent  on  the  degree  of  the  spherical  harmonic  trial 
function  expansion  which  is  used.  They  are  not  dependent 
on  the  problem  parameters  and  therefore  they  can  be  inde¬ 
pendently  evaluated.  They  can  be  evaluated  once,  and 
thereafter,  used  as  a  part  of  the  problem  input  data.  The 
source  integrals  are  derived  from  an  interpolation  of  the 
first  scatter  air-over-ground  source  over  the  entire 
spatial  problem  domain.  This  first  scatter  source  is 
derived  in  the  next  section. 
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The  Source  Terms 

A  numerical  solution  of  the  second  order  transport 
equation  requires  that  the  source  terms  (right  hand  side 
of  Eq  (126)  must  be  evaluated.  These  terms  form  the 
individual  elements  of  the  problem  source  vector  in  Eq 
(114).  The  even  and  odd  parity  sources  Q®  and  Qu  are 
comprised  of  the  relevant  terms  of  Eqs  (111)  and  (112). 

For  this  steady-state  monoenergetic  solution  these  are  the 
even  and  odd  parity  inhomogenous  sources  S®  and  Su  where 
we  have  dropped  the  explicit  multigroup  index.  These 
sources  will  be  defined  as  the  first  scatter  or  collision 
even  and  odd  parity  sources.  For  the  air -over -ground 
problem  the  first  scatter  source  S(r,fi)  is  the  number  den¬ 
sity  of  particles  which  leave  the  source  point  and  undergo 
only  one  collision  before  being  scattered  into  direction 
at  position  r.  Streaming  neutrons  which  leave  the  source 
point  and  do  not  collide  before  reaching  position  (r ,ft)  are 
not  included  in  the  collision  source. 

The  use  of  a  first  scatter  source  makes  the  air-over- 
ground  problem  more  isotropic.  It  removes  the  strongly  ani 
sotropic  streaming  particles  from  being  a  part  of  the 
problem.  Therefore,  the  solution  fluence  of  Eq  (126) 

will  be  the  scattered  even  parity  fluence  4'f(r,n)  and 
not  the  total  even  parity  fluence  ^(£,6).  Following 


the  derivation  of  Wills  (1981),  the  total  even  parity 
fluence  can  be  defined  as 


=  Vs(r,n)  +  ¥d(r,n)  (133 


where 


^.(r.fi)  =  streaming  uncollided 
particles  at  position 
(*.«>. 

A  precise  mathematical  definition  of  the  Sg  and  Su 
sources  will  now  be  developed.  Also  a  source  interpolation 
procedure  will  be  outlined.  This  source  interpolation  is 
used  in  order  to  simplify  the  source  integrals  of  Appendix 
F. 

The  First  Scatter  Source.  The  even  and  odd  parity 
sources  have  been  defined  as 


Su(r,f2)  =  %{s(r,6)  -  S(r,-6)| 


(134: 


Sg(r,6)  =  %{s(r,6)  +  S(r,-6)1 
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If  Su  and  sS  are  first  scatter  source  densities  then 
S(r,ft)  and  S(r,-fi)  must  also  be  defined  as  first  scatter 
source  particles/unit  volume.  If  a  position  (p,z)  in  the 
problem  domain  is  chosen  then  a  unit  vector  from  the 
source  point  (0,zb)  can  be  defined  as 


ft' (P ,z) 


peQ  +  (z-zb)e? 
{p2  +  (z-zb)2}% 


(136) 


Figure  8  shows  the  direction  vectors  of  this  first  scatter 
(collision)  source.  0.'  is  the  direction  that  all  streaming 
(uncollided)  particles  have  at  point  (p,z). 

By  definition  only  particles  which  are  streaming 
radially  outward  from  the  source  point  can  be  included  in 
the  direction  fluence.  Therefore  the  direct  fluence  at 
point  (p,z)  is  in  the  direction  and  can  be  written  as 

<t>d(p,z,n’)  =-5^2  exp  {-  j  ot(z)dsj  (137) 


where 


s  =  jp2  +  (z-zb)2j^ 


(138) 
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Figure  8.  First  Scatter  (Collision) 
Source  Direction  Vectors. 


and  j ds  means  that  the  integration  is  carried  out  along 
the  path  s  (See  Figure  8). 

Y  =  Total  number  of  source  particles. 

a  a  ( 0)e  z/s^  for  z  >  0 

t(z)  = 

ot  (ground)  for  z  <  0  (139) 

sh  =  atmospheric  scale  height 

The  term  ft  ot(z)ds  is  the  average  number  of  collisions 
which  a  particle  undergoes  in  traveling  from  the  source 
point  (o,  zb)  to  point  (p,z).  From  Figure  8  the  distance 
s  can  also  be  written  as 

s  =  (z-zb)/yd  (140) 

and  therefore  by  changing  variables 

ds  =  dz/yd  (141) 

where  yd  is  a  function  of  p  and  z  (but  constant  along  a 
path  length)  and, 


yd(p,z)  =  cos  (w)  =  (z-zb)/s 


(2-zb) 

]p2  +  (z-zb)2|^ 


The  integral  term  of  Eq  (137)  can  now  be  written  for 
as 


=f0S  %(»<U  -  -^£*zlshiz 


and  finally  as 

_/-  _ \  _  °t^  j  _-zb/sh  _-z/sh  | 
T(p,z)  =  --  -j-  <  e  -  e  >*sh 

-  {<*.<*»  - 

From  the  above  derivation  it  follows  that 


t(p,z) 


|at(zb)  for  z  >  0 

{°t(*b)  -  ot(0)}  ^-^forxO 


(142) 


>  0 


(143) 


(144) 


(145) 
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also 


<Mp 


.  z  .  ft  ' )  - 


V 

-  i - oext>(- 


(146) 


Note  that  <t>j(p  ,z  is  only  a  function  of  p  and  z. 

The  first  scatter  source  at  (p,z)  and  with  direction 
yd  are  those  particles  which  undergo  their  first  collision 
at  (p,z)  and  are  scattered  from  direction  to  6. 
Therefore  the  first  scatter  source  can  now  be  defined  as 


S(p,z,fi)  =  as(z,fi>n')())(j(p,z,fi')  (147) 


where  as  is  not  a  function  of  ft'  but  of  the  scattering  angle 
fi.fl'(y  )  and  z.  From  Figure  8  and  Figure  4  ft'  is  defined 
by  yd  and  x  =  0  i.e.,  ft'  =  (y',x^)  where  y'  =  yd  and 
x'  =  0. 

By  use  of  the  addition  theorem  it  is  shown  in  Appen¬ 
dix  H  that  Eq  (147)  can  be  written  as 

L  2. 

S(p,z,fi)  =  2*dG,z,b')e~z/sh'£t  2°£(0)CLPJto(yd)cos(inx)  <U8) 

2=0  m*=0 
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and  that  the  even  and  odd  parity  first  scatter  sources  are 


Sg(P,z,6) 


4>d(p  ,z,6' )  e 


+  C-D  Jl}pjlin<y)pJlTI1(>id)cos(inx) 


(149) 


and 


L  l 

Su(p,z,fi)  =  <f>d  ( p ,  z ,  ^  *  )e  ao(0)C2  /l 

£=0  m*=0  1  1  m  ( 

'  <~1)J'}pjlm(lJ)I>£ni<^)cos(inx)  (150) 

where  m*  means  that  all  terms  with  an  m  =  0  subscript  must 
be  divided  by  two. 

Source  Interpolation.  Because  of  the  complicated 
nature  of  the  source  expressions,  Eqs  (149)  and  (150),  and 
the  need  to  integrate  the  source  terms  of  Appendix  F,  a 
spatial  source  interpolation  will  be  used.  This  interpola¬ 
tion  which  simplifies  the  source  integrals  is  necessary  if 
a  very  tedious  (double  or  quadruple)  integration  is  to  be 
avoided.  By  this  interpolation  process  the  source  terms 
of  Appendix  F  can  all  be  separated  into  a  product  of  single 
integrals. 
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It  is  important  to  note  that  ud  which  is  given  by 
Eq  (142)  is  a  function  of  p  and  z.  Furthermore,  ^(p.z,^') 
of  Eq  (146)  is  a  function  of  p  and  z.  Beginning  with 
Eqs  (149)  and  (150)  they  can  be  rewritten  as 


L  £ 

Sg(p,z,fi)  Z(a®(0)<4i{l+  (-l)£}A£m(p,z)P£m(u)cos(mx) 


(151) 


and 


L  £ 


iu(P,z,S).g)^(0)c2m{l-  (-l/^p.zS^^cosCmx) 


(152) 


where 


A£m(p’z)  =  4>d(P,z,^,)P£tn(pd)e"z/sh  (153) 

A  spatial  (p,z)  interpolation  of  the  even  and  odd 
parity  sources,  Eqs  (151)  and  (152),  is  therefore  an  inter- 
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polation  of  A^m(p,z).  In  this  project  these  first  scatter 
second  order  sources  will  be  interpolated  by  a  combination 
of  piecewise  bilinear  Lagrange  polynomial  functions. 
Specifically,  S8  is  approximated  by  a  tensor  product  of 
linear  Lagrange  polynomials  as  follows 

NZ  NR 

Sg(p,z,fi)  =  2  S  S8(p.  ,z.  ,fi)H.  (p)H.  (z)  (154) 

i=l  j=i  3  1  3  1 

where 

S8(p.,zif6)  =  the  even  parity  source,  Eq  (149) 

J  evaluated  at  the  spatial  nodes 

(pj  >zi) 

NZ  =  total  number  of  z-nodes 
NR  =  total  number  of  R-nodes 
H(p)  =  p-linear  Lagrange  polynomial 
H(z)  =  z-linear  Lagrange  polynomial 

These  linear  polynomials  (See  Figure  9)  are  defined  as 


H.(x) 


xi-l  ~  x  for  x.  -i^x^x. 
xi-l  "  xi 

xi+l  ~  x  for  x.lxsx. ,, 
xi+l  ‘  xi 


(155) 
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The  product  Hj(p)-Hi(z)  of  Eq  (154)  forms  a  tensor  product 
space  on  a  rectangular  grid,  in  the  p,  z  plane  (Prenter, 
1975) 

Substituting  Eq  (149)  for  Sg(Pj,z^,6)  in  Eq  (154) 

gives 


Sg(p,z,6)  -T!  I  feojclnl1  +  (-!>*}  P^JcosOnx)* 


NZ  NR  — 

II  \m(Pj’2i)HJ(P)Hi(z)J  (156) 

i=l  j=l 


Similarly  the  odd  parity  source  can  also  be  expressed  as 


L  Z 


Su(p,z,a)  Z  k(0)C^m{l  -  (-D1Km<M)cos(mx)* 

Z=0  m*=0 


NZ  NR 


Z  Z  Ato(pj-zi)Hj<p)Hi(z)] 

i-lj-1 


(157) 


Eqs  (156)  and  (157)  can  now  be  substituted  into  the  source 
terms  of  Eq  (126).  This  substitution  and  a  separation  of 
the  integration  variables  produced  a  number  of  source  inter- 
grals.  These  integrals  are  listed  in  Appendix  G. 
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Computer  Solution  and  Results 

A  finite  element  projected  solution  of  the  air-over 
ground  problem  can  be  obtained  from  a  computer  solution  to 
a  system  of  coupled  algebraic  equations,  which  are  in  the 
form  of  Eq  (126).  In  order  to  validate,  evaluate  and 
demonstrate  the  potential  of  our  finite  element  formalism, 
a  computer  program  was  developed  to  solve  the  steady- 
state  monoenergetic  problem.  This  solution  involves  a 
computer  assemblage  and  solution  of  the  matrix  problem  of 
Eq  (126). 

Our  solution  is  basically  a  computer  implementation 
of  the  problem  equations  and  expansions  of  the  previous 
section.  We  developed  a  computer  code  to  solve  the  even 
parity  Galerkin  equation  with  vacuum  boundary  conditions, 
first  scatter  sources,  and  an  exponentially  varying  at¬ 
mosphere.  This  is  a  solution  to  Eq  (126)  with  the  source 
terms  of  Eqs  (156)  and  (157),  and  in  an  air-ground  cylin¬ 
drical  geometry  where  the  point  sources  are  confined  to  the 
z-axis.  A  tensor  product  of  cubic  splines  (k  =  A)  is 
used  to  represent  the  spatial  (p,z)  variables. 

This  is  a  four-dimensional  problem  with  two  spatial 
and  two  angular  (y,x)  variables.  The  ground  is  modeled  as 
a  homogenous  cylinder  with  constant  material  properties 
(density  and  cross-sections).  However,  in  keeping  with 
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the  exponentially  varying  density  of  the  atmosphere 
(Messier,  1971)  the  cross-section  of  air  is  modeled 
with  an  exponential  dependence  where 

cr(r)  =  a(0)e_z/H  (158 

air  air 

and 

o(0)  =  cross-section  of  air  at  sea  level 
air 

H  -  atmospheric  scale  height 
z  *  height  above  sea  level. 

In  order  to  compare  and  validate  our  solution  with  other 
methods  (codes)  we  also  included  an  option  (H  °°) 
which  models  a  homogenous  atmosphere. 

The  Computer  Program.  This  is  a  straightforward 
computer  implementation  of  the  problem  equations  which 
were  developed  in  the  previous  section.  These  equations 
are  expanded  in  Appendix  F.  The  integrals  which  must 
be  evaluated  are  in  Appendix  G. 

This  finite  element  transport  (FET)  code  computes 
the  even  parity  flux  and  the  scalar  or  total  reaction  rate 
flux  for  the  air-over-ground  problem.  It  was  developed 


specifically  to  solve  this  problem.  The  structure,  format 
and  design  of  the  FET  program  are  based  upon  the  special 
features  and  difficulties  of  the  air-over-ground  problem. 
Therefore,  the  FET  code  is  designed  for  a  specific 
research  application  and  not  as  a  general  particle  transport 
code. 

The  FET  program  consists  of  over  1800  lines  and 
twenty  subprograms .  We  used  a  modular  programming  style 
and  separated  the  computational  tasks  into  five  main  areas. 
In  each  of  these  areas  we  developed  function  and 
subroutine  subprograms  to  carry  out  specific  tasks.  We 
also  used  an  eight  point  Gauss -Legendre  adaptive  quadrature 
routine  to  numerically  evaluate  the  integrals  of  Appendix 
G.  This  integration  program  is  a  part  of  a  system 
library  of  mathematical  subroutines  and  functions.  It 
belongs  to  the  SLATEC  (Sandia-Los  Alamos-Air  Force  Weapons 
Laboratory  Technical  Exchange  Committee)  common  mathemati¬ 
cal  subprogram  library  (Allen  and  Funk,  1981). 

The  main  computational  tasks  in  the  order  in  which 
they  were  performed  by  FET  are  as  follows: 

1 .  Program  Initialization  and  Input 
The  problem  parameters  and  data  are  read  from  an 
input  file  by  the  main  program  driver.  This  pro¬ 
gram  determines  the  problem  size  and  computes  the 
problem  mesh  (elements). 
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2.  Evaluation  of  the  Integrals 

The  integrals  of  Appendix  C  are  evaluated  for  each 
of  the  trial  function  expansion  subscripts. 

3 .  Computation  of  the  Source  Vector 

The  first  collision  source  is  computed  and  the 
source  vector  (right  hand  side  of  Eq  (126))  is 
assembled.  The  individual  elements  of  the  source 
vector  are  computed  and  stored  as  a  singly  sub¬ 
scripted  real  array. 

4 .  Compute  the  Problem  Matrix 

The  individual  terms  of  the  left  hand  side  of  Eq 
(126)  are  computed.  TheF*  terms  are  then  added  and 
subtracted  to  form  the  problem  matrix.  This  is  a 
doubly  subscripted  N  x  N  array.  However,  because 
the  problem  matrix  is  positive  definite  symmetric 
only  those  elements  on  or  above  the  diagonal  are 
computed.  Furthermore,  a  special  sparse  symmetric 
storage  scheme  is  used  (Bathe,  1982),  and 
the  problem  matrix  is  stored  as  a  singly  subscripted 
array . 

5 .  Solution  of  the  Matrix  Problem 

Having  assembled  the  problem  matrix  and  source  vector 
the  FET  code  solves  a  matrix  problem  for  the  even 
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parity  expansion  coefficients  (the  A^'s)  of  Eq  (121). 
Because  the  problem  matrix  is  positive  definite 
symmetric  we  could  use  an  indirect  (iterative)  solver 
(Hageman  and  Young,  1980).  We  can  also  use  a 
Gauss -Jordan  elimination  process  or  some  other  direct 
method  (Atkinson,  1978).  FET  uses  Cholesky  decom¬ 
position  which  is  a  variant  of  Gauss-Jordan 
elimination  for  positive  definite  symmetric  matrices. 
A  subprogram  which  is  based  on  the  routine  Colsol 
of  Bathe  (1982)  is  used  in  the  FET  code.  This 
subprogram  carries  out  a  Cholesky  decomposition  of 
the  problem  matrix  and  a  back-substitution  of  the 
source  vector  (Jennings,  1977). 

Results  and  Comparisons.  The  evaluation  and  vali¬ 
dation  of  our  finite  element  projection  technique  was 
accomplished  by  using  the  FET  code  to  solve  a  number  of 
sample  problems.  These  solutions  show  that  our  solution 
technique  is  capable  of  solving  difficult  problems  of 
physical  and  engineering  interest.  They  are  not  meant  to 
be  a  complete  and  exact  solution  of  the  air -over- ground 
problem  but,  rather  a  demonstration  of  the  potential  of 
our  solution  approach  to  solve  general  particle  transport 
problems.  However,  they  do  show  that  for  the  air-over¬ 
ground  problem  our  method  can  serve  as  a  feasible  and  com- 
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plementary  alternative  to  Monte  Carlo  and  discrete 
ordinates . 

The  air-over- ground  problem  was  chosen  because  of 
its  inherent  difficulties  and  the  limitations  of  presently 
available  solution  techniques  (Loewe,  et  al.  ,  1983).  In 
order  to  validate  our  solution  approach  and  the  FET  code 
a  number  of  homogenous  air  calculations  were  performed. 
These  calculations  afforded  a  comparison  of  our  results 
to  those  of  other  methods .  We  compared  our  results  with 
those  of  diffusion  theory  and  discrete  ordinates.  We 
also  carried  out  homogenous  and  inhomogenous  P^  air- 
over-ground  calculations.  These  calculations  demonstrate 
the  potential  of  our  method  for  providing  detailed  and 
accurate  solutions  to  the  air-over-ground  problem.  We 
will  now  present  these  results  and  comparisons . 

One  group  neutron  cross-sections  from  the  DLC-31  cross- 
section  set  of  Oak  Ridge  National  Laboratory  (ORNL)  were 
used.  These  cross-sections  are  a  part  of  the  Radiation 
Shielding  Information  Center  (RSIC)  data  library.  Two 
sets  of  one-group  cross-sections  were  used  to  perform 
the  calculations.  These  cross-sections  were  assembled 
by  ORNL  from  the  Evaluated  Nuclear  Data  File  IV  (ENDF-IV). 

A  uniform  1/E  weighting  function  was  used  to  average  and 
collapse  the  cross-sections  to  37  groups. 


The  cross-sections  are  approximated  by  a  4th  order 


(P^)  Legendre  polynomial  expansion.  The  two  sets  of 
one  group  cross-sections  which  we  will  label  group  one 
and  two  reflects  the  anisotropic  (forward  peak)  scattering 
of  air.  These  cross-sections  are  for  homogenous  air  at  an 
air  density  of  1.11  mg/cc  and  with  mean  free  paths  of 
0.16  km  and  0.11  km  for  groups  one  and  two  respectively. 
The  group  two  scattering  cross-sections  are  more 
anisotropic  and  forward  peaked  than  those  of  group  one. 

In  order  to  validate  our  method  and  debug  the 
computer  code  FET  we  compared  our  results  to  those  of  an 
analytic  diffusion  theory  solution.  This  analytic  solu¬ 
tion  is  based  on  a  point  source  of  one  particle  in  a 
finite  homogenous  sphere.  It  can  be  shown  (Glasstone  and 
Edlund,  1952)  that  the  diffusion  solution  of  this  problem 
with  vacuum  boundary  conditions  is  given  by 


<Kr) 


So[e"r/L-e(r'2R°)/L] 

4ttD(1  -  e"^Ro^)r 


(159) 


where 

r  =  distance  from  the  point  source  (center  of  sphere) 
Rq  =  radius  of  sphere 
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L  =  diffusion  length 

D  =  diffusion  coefficient 
SD  =  source  strength 

Diffusion  theory  is  inappropriate  to  model  this  pro¬ 
blem  exactly,  especially  close  to  the  boundaries  and  point 
sources  (Duderstadt  and  Hamilton,  1976).  Furthermore  be¬ 
cause  the  FET  code  provides  the  solution  to  a  point  source 
on  the  z-axis  in  a  cylindrical  geometry  an  exact  comparison 
cannot  be  made.  However,  we  can  compare  these  two  results 
in  an  approximate  way. 

Both  the  FET  code  and  analytic  diffusion  results  are 
for  a  point  source  of  one  particle  and  the  group  one  and  two 
homogenous  air  cross-sections.  These  results  for  a  sphere 
and  cylinder  with  vacuum  boundaries  are  shown  in  Figures 
10  and  11.  The  diffusion  calculations  are  performed  in  a 
sphere  with  a  radius  of  2km  and  with  the  point  sources  at 
the  origin.  The  FET  calculations  are  for  a  cylinder  of 
radius  and  height  equal  to  2km  respectively.  The  source 
is  located  at  the  center  of  the  cylinder.  For  the  group  one 
calculation  the  cylinder  is  subdivided  into  256  rectangular 
elements  (cells).  However,  for  the  group  two  calculation 
it  is  partitioned  into  324  (rectangular)  elements.  These 
(FET)  results  are  presented  for  radial  points  at  the  source 
altitude. 
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To  further  examine  and  validate  our  results  we  also 
made  two  comparisons  with  discrete  ordinate  solutions. 

These  solutions  were  independently  obtained  by  using  the 
discrete  ordinates  (Sn)  code  ANISN  (Burgio,  1975).  Infinite 
homogenous  air  calculations  in  spherical  geometry  were  per¬ 
formed  with  point  sources  at  the  origin  and  an  S^2 
angular  quadrature.  These  are  published  results  for  a 
number  of  one  group  problems.  We  have  presented,  in 
Figures  12  and  13,  a  comparison  of  these  results  for  the 
cross-sections  of  groups  one  and  two.  The  FET  results 
are  for  a  2km  x  2km  homogenous  cylinder  with  vacuum 
boundary  conditions.  Here  again  the  cylinder  is  parti¬ 
tioned  into  256  and  324  rectangular  elements  for  the  groups 
one  and  two  calculations  respectively.  The  point  source  is 
at  the  cylinder  center  and  the  flux  values  are  for  radial 
points  which  are  at  the  source  altitude.  Because  of 
the  differences  in  the  geometry  of  the  FET  and  Sn  cal¬ 
culations  an  exact  comparison  cannot  be  made.  However,  we 
are  able  to  provide  an  approximate  comparison  and  thus  demon¬ 
strate  that  our  results  are  appropriate  and  valid. 

We  also  performed  a  number  of  air- over -ground 
calculations.  The  homogenous  ground  cross-sections  for 
group  one  and  two  are  larger  and  less  anisotropic  than 
for  that  of  air.  Furthermore,  the  ground  cross-sections 
for  group  two  are  more  anisotropic  than  the  group  one 
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Figure  13.  Group  2  --  Discrete  Ordinates  Comparison 


ground  cross-sections.  The  ground  mean  free  paths  for 
group  one  and  two  are  1.05E-4  and  7. 42E-5  kilometers .  At 
the  air-ground  interface  the  material  properties  and 
cross-section  data  are  rapidly  varying.  The  air  and 
ground  cross-sections  vary  by  three  orders  of  magnitude. 

This  presents  special  difficulties  and  a  severe  test  for 
most  numerical  solution  techniques  (Lowe,  et  al. ,  1983). 

The  difficulties  in  solving  the  air-over-ground 
problem  are  also  due  to  the  presence  of  point  sources  and 
an  exponentially  varying  atmosphere.  These  difficulties 
are  further  compounded  by  the  need  to  properly  treat  the 
air-ground  interface  condition.  This  is  an  important 
consideration,  especially  when  a  solution  is  required  at 
the  interface.  For  many  applications  a  solution  is  required 
at  ground  level  or  very  close  to  the  interface  (Pace,  et 
al. ,  1975).  Therefore,  we  are  primarily  interested  in  a 
solution  to  the  air-over-ground  problem  for  radial  points 
(ground  range)  along  this  interface. 

This  interface  condition  is  treated  as  a  natural 
condition  in  our  finite  element  projected  solution.  There¬ 
fore,  our  solution  also  demonstrates  the  ability  of  our 
method  to  model  the  air-ground  interface.  Solutions  for 
the  homogenous  and  exponentially  varying  problems  are 
presented  in  Figures  14  through  17.  The  group  one  and  two 
cross-sections  were  used  to  calculate  scalar  flux  values 
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Figure  15.  Group  2  — Homogenous  Air-Over-Ground  Solution 
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Figure  16.  Group  1  --Exponential  Air-Over-Ground  Solution 
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Figure  17.  Group  2 


Exponential  Air -Over -Ground  Solution 
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at  the  air-ground  interface.  The  problem  domain  is  a 
2km  x  2km  cylinder  with  vacuum  boundary  conditions.  The 
cylinder  has  304  cells  for  the  group  one  problem  and  381 
cells  for  group  two.  A  scale  height  (H)  of  7km  is  used 
in  the  exponentially  varying  air-over-ground  calculations. 

In  Figures  18  and  19  we  show  the  variation  of  the  group  one 
homogenous  air-over-ground  solution  across  the  interface. 
This  (air-ground)  interface  is  located  at  4.5E-4  kilometers 
from  the  bottom  of  the  problem  cylinder. 

Finally,  we  examined  contour  plots  for  the  homo¬ 
genous  air  (only)  problem  and  the  convergence  of  the  FET 
solution.  Figure  20  shows  an  example  of  this  convergence. 
These  are  results  at  points  along  the  air-ground  interface 
for  the  group  one  homogenous  air -over -ground  problem.  The 
problem  domain  is  a  1km  x  1km  cylinder  which  is  partitioned 
into  a  nonuniform  mesh  with  an  equal  number  of  mesh  points 
along  the  z  and  p  axes. 

In  Figure  21  we  present  a  contour  plot  of  the  group 
one  FET  homogenous  air  solution.  This  plot  is  intended 
to  show  any  deviation  of  the  solution  from  the  physical 
symmetries  of  the  problem.  Such  deviations  or  distor¬ 
tions  of  the  solution  is  usually  an  indication  of  ray 
effect  in  the  discrete  ordinate  solution  (Loewe,  et  al., 
1983).  Figure  21  shows  that  the  FET  solution  does  not  suffer 
from  the  anomalous  ray  effect  problem  of  discrete  ordinates. 
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These  steady-state  one  group  results  show  the  poten¬ 
tial  of  this  solution  technique  to  provide  a  complete 
solution  to  the  air-over-ground  problem.  However,  an 
energy  dependent  multigroup  solution  is  required  for  most 
applications.  For  others  a  complete  time  and  energy 
dependent  solution  will  be  needed.  In  either  case,  these 
solutions  can  be  obtained  by  solving  a  number  of  coupled 
one  group  problems.  Therefore,  our  steady-state  one  group 
solution  can  be  extended  to  the  general  air-over-ground 
problem. 

All  of  our  computations  were  carried  out  on  a  Cray- 
I  computer.  The  FET  code  was  written  in  Fortran  77  and  it 
is  not  optimized  in  terms  of  storage  requirements  or 
efficiency.  The  calculations  were  all  carried  out  in- 
core  (central  memory) .  The  computer  storage  requirements 
were  such  that  special  auxiliary  or  disc  storage  was  not 
needed.  Furthermore,  the  maximum  execution  time  which 
was  required  was  less  than  ten  minutes.  However,  for 
larger  problems  we  expect  the  storage  requirements  and 
execution  times  to  increase. 


CHAPTER  IX 


CONCLUSION  AND  RECOMMENDATION 

A  method  for  solving  general  particle  transport 
problems  has  been  developed.  This  method  is  based  on  a 
mathematical  formalism  which  includes  a  rigorous  treatment 
of  general  boundary  conditions .  These  boundary  conditions 
are  dealt  with  within  the  context  of  a  classical  energy 
minimization  (variational)  principle  and  an  equivalent 
Bubnov-Galerkin  solution.  Using  the  second  order  forms 
of  the  transport  equation,  this  equivalence  is  established 
and  the  resulting  matrix  problem  is  positive  definite 
symmetric . 

The  mathematical  formalism  also  includes  an  explicit 
forward  or  weighted  difference  in  time  and  a  piecewise 
constant  energy  dependence.  This  is  the  usual  multigroup 
treatment  of  the  energy  dependence.  With  these  approxi¬ 
mations  of  the  time  and  energy  dependences,  the  transport 
problem  (equations)  is  reduced  to  a  set  of  coupled  steady- 
state  monoenergetic  problems. 

The  numerical  technique  involves  a  finite  element 
projected  solution  of  the  transport  equation.  A  finite 


element  trial  space  of  spherical  harmonics  and  polynomial 
splines  is  used.  The  particle  flux  is  expressed  as  a 
linear  and  separable  sum  of  even  and  odd  components  on  the 
unit  sphere.  Then,  a  numerical  solution  using  the  parity 
(odd,  even)  characteristics  of  spherical  harmonics  was 
developed.  In  this  development  we  included  anisotropic 
sources  and  scattering.  Furthermore,  we  established  that 
in  the  diffusion  limit  our  method  reduces  to  the  method 
or  diffusion  theory. 

A  validation  of  the  method  has  been  obtained  by  a 
computer  solution  to  the  air -over -ground  problem.  This 
problem  is  modeled  in  cylindrical  (r,z)  geometry  with  an 
exponentially  varying  atmosphere,  anisotropic  scattering 
and  anisotropic  first  scatter  sources.  The  ground  was 
treated  as  a  homogenous  mixture  and  vacuum  boundary 
conditions  were  used.  The  computer  code  FET  was  developed 
and  used  to  solve  the  problem. 

Comparisons  of  the  FET  solutions  to  those  of  dif¬ 
fusion  theory  and  discrete  ordinates  validated  our  solution 
technique.  These  comparisons  show  the  ability  of  our 
method  to  solve  general  and  difficult  particle  transport 
problems.  Furthermore,  the  potential  of  this  solution 
approach  to  serve  as  a  complementary  alternative  to  the 
standard  techniques  of  Monte  Carlo,  discrete  ordinates  and 
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the  Pn  method  was  demonstrated.  However,  the  use  and 
validation  of  this  method  by  solving  other  particle  trans¬ 
port  problems  is  required. 

Further  work  is  also  required  in  the  computer 
implementation  of  this  method.  The  applications  of  this 
technique  to  general  time  and  energy  dependent  problems  is 
required.  A  number  of  different  finite  element  types  and 
shapes  should  be  used  and  examined  as  to  their  potential 
to  model  other  problems  of  physical  and  engineering  interest. 

Because  of  the  special  properties  of  the  problem 
matrices  the  method  should  be  examined  in  terms  of 
computational  efficiency  and  costs.  Furthermore,  other 
comparisons  should  be  made.  Comparisons  for  problems 
where  solutions  by  other  methods  are  available  and  can 
be  easily  obtained  are  required.  For  problems  with  angular 
singularities,  the  issues  of  convergence  and  stability 
of  the  solution  should  be  examined. 

In  this  regard,  we  feel  that  our  method  can  be 
easily  applied  to  problems  with  difficult  and  complicated 
geometries.  Nuclear  reactors  and  radiation  shielding 
problems  are  in  this  category.  However,  many  of  these 
problems  have  singular  points.  These  are  points  where  the 
solution  has  unbounded  derivatives  (Strang  and  Fix,  1973). 
Furthermore,  these  angular  singularities  may  reduce  the 
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accuracy  and  rate  of  convergence  of  our  finite  element 
method. 

The  equivalence  of  this  technique  to  diffusion  theory 
should  be  further  examined.  In  the  limit  of  diffusion  the¬ 
ory,  this  method,  with  a  finite  element  basis  of  linear 
splines,  should  produce  the  diffusion  theory  finite  dif¬ 
ference  matrix.  This  is  a  Stieltjes  matrix  which  guarantees 
an  all  positive  solution.  Furthermore,  a  solution  by  other 
projection  methods  such  as  collocation,  least  squares  and 
the  method  of  moments  should  be  examined.  Collocation  is 
especially  attractive  in  that  there  exists  an  equivalence 
between  the  Galerkin  and  collocation  methods  (Prenter, 

1975). 

Finally,  a  Double  Pn  aproximation  should  be  applied 
to  the  air-over-ground  problem.  The  effects  of  this 
approximation  on  the  solution  accuracy,  especially  at  the 
air-ground  interface,  should  be  examined.  This  approxi¬ 
mation  would  also  be  able  to  better  model  the  vacuum 
boundary  conditions.  However,  the  determination  of  a 
local  reflection  coefficient  and  use  of  the  albedo  boundary 
condition  might  be  more  appropriate.  Furthermore,  the 
applicability  of  our  solution  technique  to  specific  charge 
particle  transport  problems  should  be  investigated.  This 
would  entail  the  inclusion  of  the  Fokker-Planck  collision 
terms  into  the  mathematical  and  numerical  models. 


APPENDIX  A 


VARIATIONS  OF  A  FIRST  ORDER  FUNCTIONAL 


First  order  functionals  provide  the  variational  model 
for  many  physical  systems.  These  functionals  which  are 
usually  quadratic  (as  in  Eq  (53)),  can  be  used  to  provide 
solutions  to  a  large  number  of  problems.  These  solutions 
are  obtained  by  determining  the  stationary  point  or  extre¬ 
mum  of  the  functional.  The  classical  Raleigh-Ritz  method 
is  usually  used  to  numerically  determine  this  extremum. 

However,  before  obtaining  a  numerical  solution  we 
are  oftentimes  interested  in  determining  the  analytic  con¬ 
ditions  that  would  make  the  functional  stationary.  These 
conditions  are  obtained  by  using  the  calculus  of  variations 
to  extremize  the  functional.  We  are  required  to  take 
variations  (derivatives)  of  the  functional  which  in  turn 
provide  us  with  information  about  the  problem  and  its  solu¬ 
tion. 

Following  the  discussion  of  Chapter  V  we  will  begin 
with  a  first  order  functional  (Gelfand  and  Fomin,  1963). 

Ku)  =  JTf(  r,u(r),Vu(r))dr  (160) 


where  V  is  the  usual  gradient  operator  in  some  n- dimen¬ 
sional  phase  space  and  we  note  that  the  variational  problem 
is  independent  (invariant)  of  the  coordinate  system  (Cou- 
rant  and  Hilbert,  1953). 

Based  on  our  earlier  definition  of  a  variation,  we 
can  write  the  first  variation  of  Eq  (162)  as 

3I(u)  =  /  3F(r,u(r) ,V  u (r))  dr  (161) 

J  R 


and 


3F  -  F(r,u*(r) ,Vu*(r))  -  F(r ,u(r) , Vu(r) )  (162) 


where 


u*(r)  =  u(r)  +  £({)(r) 


(163) 


and 


3u(r)  =  e<j)(r) 


(164) 
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Expanding  9F  in  a  Taylor  series  and  keeping  terms  of 
order  e;  or  taking  the  variational  derivative  of  F  (Lancoz, 
1949)  we  get 


3I(u) 


(165) 


where  we  have  dropped  the  explicit  r  dependence  and 
u'  *  Vu  ;  <f>'  =  V 4> . 

For  a  stationary  point  we  require  that  Eq  (165)  be 
equal  to  zero.  Then  we  can  expand  the  individual  terms  of 
the  integral  and  discuss  the  conditions  for  an  extremum  to 
exist.  However,  to  decide  on  whether  Eq  (160)  is  a  minimum 
or  maximum  principle  we  must  take  the  second  variation. 

This  is 


32I(u)  =  3(3I(u)) 


(166) 
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APPENDIX  B 


THE  EULER-LAGRANGE  EQUATION 

The  Sturm-Liouville  variational  problem  of  Chapter  V 
corresponds  to  a  linear  boundary  value  problem  with  mixed 
boundary  conditions.  This  is  a  natural  boundary  condition 
of  the  variational  or  equivalent  Galerkin  problem.  In 
order  to  show  the  equivalence  between  the  variational  and 
boundary  value  problems  we  will  take  the  first  variation  of 
Eq  (53).  This  procedure  will  produce  the  Sturm-Liouville 
equation  and  identify  the  natural  boundary  conditions. 

In  the  terminology  of  variational  calculus  the  partial  dif¬ 
ferential  equation  which  extremizes  the  functional  is 
called  its  Euler-Lagrange  equation.  For  the  model  problem 
of  Chapter  V,  this  is  (the  Sturm-Liouville  equation)  Eq 
(50). 

We  begin  by  taking  the  first  variation  of  Eq  (53)  in 
accordance  with  the  derivation  of  Appendix  A.  This  gives 


9I(u) 


zj  2  p(r)Vu(r)4>’  +  £q(r)u(r)  -  f(r)J<Kr)  dr 


+  c(rs)]<J>(rg)ds 


(168) 


and  since  <j>'  =  V4>(r)  we  can  use  the  product  rule  of 
differentiation  and  the  Gauss  divergence  theorem  to  get 


3I(u)  =  ejT  2  |-V(p(r)  Vu(r))  +  q(r)u(r)  -  f(r)j<J>dr 


tf  2  fp 

•'s  L 


+  e£  2  jp(rg)Vu(rs) -n  +  b(rg)u(rg)  +  c(rg)  |(()(fg)ds 

(169) 
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where  n  is  an  outward  unit  normal  to  the  problem  boundaries. 

For  a  minimum  to  exist  we  require  that  the  first 
variation  should  be  zero  or  that  the  variational  derivative 


3I(u)  =  lim  9I(u)  (170) 

9e  e+0  e 


be  zero.  In  either  case,  since  e  is  small  and  not 
equal  to  zero  and  d> (or )  is  an  arbitrary  function,  we  must 
have  that  4>(r)  is  equal  to  zero  everywhere  or  that 


-V(p(r ) Vu(r) )  +  q(r)u(r)  -  f(r)  =  0 


(171) 


and 


p(rg)Vu(rg) *n  +  b(rg)u(rg)  +  c(rg)  =  0  (172) 


Eq  (171)  is  the  Sturm-Liouville  problem  of  Eq  (50). 

It  is  also  the  Euler-Lagrange  equation  of  Eq  (53).  Further¬ 
more,  Eq  (172)  represents  the  natural  boundary  conditions 
which  are  satisfied  automatically  by  the  minimizing 
sequence.  A  numerical  solution  of  the  Raleigh-Ritz  system 
of  Eq  (61)  will  provide  a  solution  to  Eqs  (171)  and  (172) . 
This  is  guaranteed  if  the  solution  is  also  a  minimum  of 
the  functional  of  Eq  (53).  Therefore,  we  have  shown  the 
equivalence  between  the  variational  and  boundary  value 
problems  of  Eq  (53),  (50)  and  (51). 

It  is  important  to  note  that  for  the  Dirichlet 

boundary  condition  we  set  a(r  )  to  zero  in  Eq  (51)  and 

s 

then  Eq  (172)  is  satisfied  if  p(r  )  is  also  equal  to  zero. 

o 

Since  p(r  )  must  be  nonzero  in  the  Sturm-Liouville  problem 
s 

then  we  must  require  that  u(rg)  be  a  constant  on  the 
boundary.  This  implies  that 


b(rg)u(rg)  =  -c(rg)  =  c  (173) 
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and  that 


where 


*(rs)  =  0 


c(rg)  is  a  constant  and  there  is  no  variation  of 
on  the  boundary. 


(174) 
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APPENDIX  C 


THE  EVEN  PARITY  FUNCTIONAL 

Based  on  the  self-adjoint  nature  of  the  second  order 
transport  equation,  we  can  formulate  quadratic  functionals 
whose  minimum  will  provide  solutions  to  general  particle 
transport  problems.  Such  a  functional  in  the  form  of  Eq 
(62)  has  been  proposed.  An  important  aspect  of  this  formu¬ 
lation  is  the  treatment  of  general  boundary  conditions 
as  natural  conditions.  In  order  to  show  that  the  conditions 
whereby  Eq  (62)  has  a  minimum  are  also  a  solution  to  the 
even  parity  transport  problem,  we  must  take  the  first 
variation  of  this  functional.  Following  the  arguments 
of  Appendices  A  and  B,  we  begin  by  writing  the  first 
variation  of  Eq  (62)  as 


3I(u) 


3F(r  ,fi  ,u,V u)dp 

8L 


(175) 


/V 

where  we  have  dropped  the  explicit  (r  ,&)  dependence  and 

u  =  ;  u'  *  V  V®  =  ;  <f> '  =  0*V<f>  and  the  integral 

Si  SL  & 

is  to  be  carried  out  over  the  entire  volume  in  phase  space. 
We  note  that  the  parity  operators  are  positive  definite 
and  self-adjoint,  where  for  a  self-adjoint  operator  L 


<f,Lh>  -  <Lf,h> 


(176) 


where  f  and  h  are  continuous  and  differentiable  functions. 
Then,  after  some  algebra  and  use  of  the  divergence  theorem 


f <ft*Vf,h>  =  -/  <f,n*Vh>dr  +  p*<(fi-n)f ,h>ds  (177) 
•7R  *7 R  J s 


we  get 


4»  -  2 <Q8 , <f> >  +  2<n*VKuQu,<{>> 


fi  -  [2<gV, 

-  2<6*VKu(fi-V'i'g)  ,<}»jdr  +  [<  |S*n|  (Qeb  +  a^8)  ,4>> 


-  2<(fi 


n)4'U,<J>>j 


ds 


(178) 


where  we  have  used  the  fact  that 


3Qeb 

3u 


a 


(179) 


and  earlier  defined  to  be  the  odd  parity  flux. 

That  is 


Vu(r8,fi)  -  Ku[Qu(rs,n)  -  6*  VH^r^)]  (180) 


We  can  now  rewrite  the  surface  term  of  Eq  (178)  as 


j(«8-fi)[qeb  +  afS  -  2l'u],*>fi.ft>0da 
+  /«“-">[Qeb  +  afg  +  <  0dS  (181) 


and  if  $  is  nonzero  we  must  require  that 


-6«VKu(fi-V4'8)  +  G8'!'8  -  Q8  +  6-VKuQu  *  0  (182) 
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and  also  that  Eq  (181)  be  equal  to  zero  in  order  for  a 
minimum  to  exist.  Therefore  Eq  (62)  has  a  minimum  if 
Eq  (178)  is  equal  is  zero.  The  conditions  for  this  to  occur 
are  that  the  Euler-Lagrange  equation,  Eq  (182),  and  the 
boundary  term,  Eq  (181),  should  be  zero.  Eqs  (182)  and 
(181)  represent  the  even  parity  transport  equation  and 
boundary  conditions  of  Chapter  IV. 

We  can  easily  recognize  that  the  Euler-Lagrange 
equation  is  just  Eq  (35).  Furthermore,  if  <j>  is  nonzero 
on  the  problem  boundaries,  then  from  Eq  (181)  we  must  re¬ 
quire  that 

Qeb  +  a¥g(rs,fi)  -  2fu(rg,n)  *  0  for  fl-n>0  (183) 


and 


Qeb  +  a'}'g(ra,n)  +  2YU(rg,n)  *  0  for  fi-n<0  (184) 
With  Qgb  defined  (See  Chapter  V)  as 


a1»g(fs,6)  -  bqg(rg,n) 


(185) 


where 


1  -  ot(rg) 
a  ■  + 1 

and  0<a(rg)<l.  We  can  now  recognize  that  the  natural 
boundary  conditions  are 

1 .  Vacuum  boundary 
a(rg)  =  0  5  b  -  0 

2.  Incident  source 
a(rs)  =  0  ;  b  *  2 

3 .  Albedo  condition 

0  <  a(rg)  <  1  ;  b  =  0 

and 

4.  Mixed  condition 

0  <  «<rs)  i  1  i  b  - 

s 

However,  the  Dirichlet  condition  with 
0  <  a(t,)  <  1  ,  b  ■  0  or 

s 


(186) 
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and  the  angular  flux  completely  or  partially  specified  on 
the  boundaries  is  an  essential  condition.  This  condition 
must  be  imposed  on  the  trial  space  in  order  to  insure  that 
it  will  be  included  in  a  solution  to  the  variational  problem, 
Eq  (62).  There  is  some  flexibility  in  specifying  the 
Dirichlet  condition.  Nonetheless,  they  must  be  consistent 
and  reflect  the  physics  of  the  problem. 

It  is  important  to  note  that  the  vacuum  boundary 
is  really  a  special  case  of  the  Dirichlet  boundary  condi¬ 
tion.  Furthermore,  if  we  use  a  global  basis  in  angle 
(spherical  harmonics)  we  are  unable  to  explicitly  enforce 
the  boundary  condition  as  essential.  Nonetheless,  we  are 
able  to  treat  the  vacuum  boundary  as  natural  with  a  *  b  *  0. 
This  approach  which  is  an  approximation  to  the  true  vacuum 
boundary  results  in  the  usual  Mar shack  boundary  condition 
(Davis,  1966). 

We  can  also  treat  the  usual  interface  conditions  with¬ 
in  this  variational  context  (Strang  and  Fix,  1973).  Noting 
that  the  variational  problem  is  dependent  upon  material 
properties  (cross-sections),  we  divide  the  problem  domain 
into  regions  of  constant  or  continuously  varying  cross- 
section  data.  In  each  region  we  require  that  the  variational 
problem  and  the  second  order  transport  equations  be  satis¬ 
fied.  We  then  treat  discontinuities  of  material  properties 
as  a  set  of  coupled  boundary  value  problems. 
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Except  for  those  problems  with  interface  singulari¬ 
ties  (Strang  and  Fix,  1973)  we  can  proceed  by  imposing  an 
albedo  boundary  condition  at  the  interface.  At  interfaces 

we  set  a(r_)  =  1  and  b  »  0  in  the  functional  of  Eq  (62). 
s 

We  then  write  a  separate  functional  for  each  region  and 
require  that  the  first  variation  vanishes  on  the  entire 
problem  domain.  The  boundary  conditions  are  treated  in 
the  usual  manner  as  essential  or  natural  conditions.  How¬ 
ever,  upon  taking  the  sum  of  the  individual  functionals  we 
require  that  the  interface  terms  vanish. 

The  interface  terms  can  be  written  as 

y  <  J(fl-n)¥]J(fI,fi)  +  (-n-n  )V^(rx,-fl)  J , <J>>ds  (187) 

where  +  means  in  the  limit  approaching  the  interface  in  a 
given  positive  0,  direction  from  the  left  and  -  means  in  the 
opposite  -fi  direction  from  the  right.  Then  with  h  being 
the  outward  normal  to  the  interface  we  have 

n_  »  -n+  (188) 

and 

(S-n+)  «  (-fi-n_)  (189) 
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We  can  now  rewrite  Eq  (187)  as 


<(fi*n+) 


[*+(*! »fi)  +  »*>d® 


(190) 


and  require  that  for  Eq  (189)  to  hold  for  arbitrary  <f> 


¥j(rIfft)  =  -Vu(rIf-fi) 


¥u(rj,S) 


(191) 


where  we  have  used  the  odd-even  properties  of  the  odd  parity 
flux. 

Furthermore,  upon  recognizing  that  4>  represents  the 
even  parity  flux  V®  which  is  an  even  function  on  the  unit 
sphere  and  fi*n  is  always  odd,  we  can  write  Eq  (190)  as 


/<(fi-n,)M'  (r_,6)  -  ¥  (rT  ,  6)1 ,4>>ds  =  0 


(192) 


(193) 


4'(rI,fi)  =  +  Vg(rr,fi) 

Then  we  get 


¥+(rI,fl)  =  H'_(rI,n) 


(194) 


Eq  (192)  requires  that  the  angular  flux  be  continuous 
at  the  interface.  This  is  a  continuity  requirement  in 
r  for  a  given  direction  which  also  implies  continuity  of 
the  particle  current  across  the  interface.  The  usual 
interface  condition  is  therefore  satisfied  as  a  natural 
condition  of  the  variational  problem. 


i 


186 
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APPENDIX  D 


THE  DIFFUSION  LIMIT 


In  Chapter  VII  we  presented  a  finite  element  pro¬ 
jection  of  splines  and  harmonic  functions  in  a  solution 
strategy  for  solving  the  second  order  transport  equation. 
In  this  Galerkin  solution  of  Eqs  (35)  and  (36)  it  can  be 
shown  that  in  the  limit  of  diffusion  theory  our  equations 
reduce  to,  and  are  consistent  with  the  diffusion  approxi¬ 
mation.  The  even  parity  Galerkin  equation  is  just  the 
diffusion  equation  whereas  the  odd  parity  equation  reduces 
to  Fick's  law. 

In  the  P^  (diffusion)  approximation  we  assume  that 
all  sources  are  isotropic  (Duderstadt  and  Hamilton,  1976) 
and  that  the  flux  and  scattering  are  linearly  anisotropic. 
Then  the  even  parity  flux  and  the  net  current  can  be 
written  as 


(195) 


J(r)  =  f  (fi*n)'l'u(r ,fi)dn 
M 


(196) 


! 


where,  in  rectangular  coordinates  (See  Figure  3) 


fl-ft  =  n-[ex  +  gy  +  6Z] 


■ryj- 

LmJ  Cn 

n  —1 _ n  "Ul 


£=lm=0 


i 

The  parity  sources  are 


Qg(r ,Q) 


S(r) 

~nr 


and 


Qu(r ,n)  «  0 


(197) 


(198) 


(199) 


Then  by  Eqs  (39)  and  (42)  and  with  the  usual  cross- 
section  expansion  (Wills,  1981)  we  can  carry  out  the  alge¬ 
bra  to  write  the  Galerkin  projection  in  angle  of  Eq  (35) 


£[l«9>  -  s<«(<&)i]i&dii 


where 


tR 


*  -  O-i 


=  0*.  -  a. 


and  o-^  and  oQ  are  the  Legendre  cross-section 
coefficients . 


(200) 


£)’  <201> 


(202) 


(203) 


expansion 
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Upon  evaluation  of  the  integrals  in  Eq  (200)  we 


get 


-V-DV<j>(r)  +  ar<p  (r)  =  S(r)  (204) 


where 


D 


=  diffusion  coefficient 


(205) 


We  can  also  rewrite  the  odd  parity  equation,  Eq(36)  as 
(Wills,  1981) 


11  A  11  A  A  11  A  A 

Gu(r)fu(r,n)  =  Qu(r,n) 


A  CT  .  A  A 

n*V'{'fe(r,fi) 


(206) 


and  carry  out  the  Galerkin  projection  on  this  equation  to 
get 


J  (r)  =  -DVcf>(r) 


(207) 


190 


Eq  (204)  is  the  diffusion  equation  and  Eq  (207)  is  Fick's 
law.  Therefore ,  for  a  linearly  anisotropic  approximation 
in  angle  the  Galerkin  projection  of  the  second  order 
transport  equation  is  consistent  with  or  diffusion 
theory . 
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APPENDIX  E 


SURFACE  HARMONIC  TENSOR 

The  use  of  surface  harmonic  tensors  in  solving 
particle  transport  problems  is  not  a  new  concept.  The 
PR  method  uses  spherical  harmonic  functions  to  expand 
and  evaluate  moments  of  the  first  order  transport  equation. 
The  scattering  cross-sections,  gradient  and  surface  normal 
operators  are  usually  expanded  in  spherical  harmonics. 
Therefore,  in  a  numerical  solution  of  transport  problems, 
the  properties  and  relationships  of  harmonic  functions 
play  a  vital  role. 

In  our  second  order  formulation  of  Chapter  VII  we 
have  found  it  necessary  to  develop  and  use  these  properties. 
Some  of  the  most  useful  ones  will  now  be  presented.  The 
harmonic  functions  are  as  defined  earlier  and  in  Chapter 
VII.  We  begin  by  writing  the  usual  scattering  cross- 
section  and  surface  normal  expansions  as 


£  r  i 

Os(i.fi-S')  =  2^ol(r)^  km(S')(4n(S>  +  )Qtm<fi) 

£-  0  m*=CL  J 


(208) 


and  in  rectangular  geometry 


(209) 


/v 


where  the  *  in  m*  means  to  divide  by  two  when  m  *  0  and 
we  have  used  the  addition  theorem  (Bell  and  Glasstone,  1970) 
to  derive  Eq  (208).  Using  Eq  (208)  we  can  then  rewrite  the 
second  order  operators  in  terms  of  Q^'s  (Wills,  1981). 

We  can  also  write  the  streaming  operators  in  terms 
of  surface  harmonic  tensors.  An  example  of  this  is  in 
rectangular  coordinates  where  we  can  write 


fi*Vf (r,«) 


3 

lx 


3_  +  Qlo  3_ 
9y  3z 


•  f(r,Q)  (210) 


Then  in  the  Galerkin  weak  form  we  can  use  the  following 
relationships  to  expand  and  simplify  the  equations 


Q10(^Qtoi(^  =  ITTT 


(ft-m  +  1) 


C  C  c 

£m  0c  j.  /0  to i  Q 

err  +  <Um)cr,  1 

Jt+l  ,m  £-l,m 


tol,m 


(211) 
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'8+l,m 


Wm 

2(28+1) 


1  0C 
m+1  ^8+1  ,m+l 


8+1 


(I 


8,-1  ,m 


:8-l,tn 


(21 


1  0c 

m+I  ''8.-1  ,m+l 


8-1 


m  +  1)(  8  +  m  +  2) 

-  ^8+1 ,m~l 


m-1" 


.  (8+m-l)(8+m)rc 
+  — ^ 


m- 


8+1 


'8-1 


( 


*  2(2UI5  “m+T  Qfc+l,m+l  '  m+I  QA.1>m+1 


CllC«m 


_ i 

m+ 

.c*+ 


q-m+1)  q-m+2) 


c  -  a+m-1)  q+m),jc 
^+l,tn-l  ^ , m- 1  ~1  ,m"l. 


To  compute  the  scalar  flux  and  total  current  we  can 
use  the  following  orthogonal  properties.  Let 


co  l  r 

g(6) 

£=0m*=0  L 


Q?«<S)  +  QL<«> 


then 


f  g(H)d$  =  V^" 

J'* IT 


J g(fi)Q^m(fi)  dQ  =  1.0  for  m  -  0 


A(fl)Qj1,(8>dS 


%  for  m  ^  0 


APPENDIX  F 


THE  AIR-OVER-GROUND 
PROBLEM  EQUATIONS 


In  order  to  obtain  a  numerical  solution  to  the  air- 
over-ground  problem  it  is  necessary  to  expand  the  Galerkin 
weak  form  of  the  even  parity  transport  equation.  The 
trial  and  weight  functions  of  Chapter  VII  will  be  used  in 
this  expansion.  Because  this  is  a  very  tedious  and  lengthy 
derivation,  the  more  obvious  algebraic  steps  will  be 
omitted.  The  starting  point  of  this  formulation  is  Eq 
(126)  of  Chapter  VIII. 


IZ  IR  L  £ 

IIIlA 

iz=lix=0  £=1  ofO 


ij£m|^  |<n-V(B; 


jQkn)'KU(^V(BiQ£m))> 


+  <BiQkn«Gg<Bi 


r  /* 


<6*V(BjQkn),KUSU>  + 


<BJ‘Vsg> 


(220) 


where  and  Q^m  are  defined  by  Eqs  (127),  (128)  and 

(96),  and  2,  and  k  are  restricted  to  be  even  integers. 

The  surface  normal  terms  are  defined  by  Eqs  (130) 
and  (131)  and  we  have  rewritten  the  source  terms  Q®  and  Qu 
as  S®  and  Su.  These  are  given  by  Eqs  (156)  and  (157). 

The  second  order  operators  Ku  and  K®  can  also  be  written 
as 

Kuf (H)  =  o.”1  ff (Q)  -  2^  T 

L  £=0m^0 

f  (fiT )dQ'|  (221) 


ot-o^ 


and 


Ggf (6)  «  otf(fi)  - 

£=0m*=0 


T2m<n,)f(n')dfi’ 


(222) 


and 


Tnm(S'>  -  +  Qtm<M'.x')Q“n(M.x)  (223) 


g 

and  Q&m  is  defined  by  Eq  (97).  Furthermore,  we  can  rewrite 
the  directional  derivative,  Eq  (126)  as 


-  fljCp*) 


I  9  l<j>B(y  x)  i  +  y 

p  ax  az 


(22  4) 


where 


a 


=  a(y,x)  =  y  1  -  y2  :os  x 


(225) 


and 


B  =  B(y,x) 


sin 


X 


(226) 


Using  these  relationships  and  noting  that  , 

®i»  Bj ,  etc.  are  all  real  functions  we  can  expand  Eq  (220) 
to  get 


IZIR  L  i 

III  Z\j,i1>n 

iz-lir-12^0  m=0 


f  -1  3  /  b  J,  dv/aaQ®  Q‘ 

JRat  3-p(PBi>72  Xr 


kn 


(227) 
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w 


+ /R&  hhdvlkm^\&^  <&'»* '  I dS 

-  ^fjx' Kn>|A<4>TrsdS'  |dfi 

+  >HPQknfcQ“T-dS  ’  H 

-  ii‘Ej>diu^((^(6qC^T«d6,  H 

+  $r  4<Bi>li<Bj)d^uQto{^Trsdfl'}d(i 


+  /0tBjBidx4QknQMd[! 


-  *t  t  /R^BjBid^knjfcsdSjdS 

r=0s*-0  R  m  <*”  '  ■ 

+  /BiBjds/  lugryi  -n2  cos  X  I  <5knQ«mdn 

S  |»IT 

■/R°‘t^lp(PB. j)dvX<nsUd“ 

-  f°l'  ^-ivfjk* to"8"" 

+  /«;'  fj<Bj>dv/<nSudS 


(239) 

(240) 

(241) 

(242) 

(243) 

(244) 

(245) 

(246) 

(247) 

(248) 

(249) 

(250) 


I 
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(251) 

(252) 

(253) 

(254) 


where 


(255) 


and 


dydx  (256) 

Eqs  (227)  to  (244)  is  an  expansion  of  the  first  term  of  Eq 
(220).  Eqs  (245)  and  (246)  are  expansions  of  the  second 
term.  Eq  (247)  is  an  expansion  of  the  third  term.  Eqs 
(248)  to  (254)  is  an  expansion  of  the  right  hand  side  of 
Eq  (220).  a®,  ct  and  c^r  are  functions  of  z  and  they  must 
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be  included  in  the  spatial  or  dv  integrals .  In  cylindrical 
geometry  with  azimuthal  symmetry 

dv  =  2irrdrdz  (257) 


and  ds  means  an  integration  over  the  surface  of  the  problem 
cylinder,  Figure  7. 

For  the  air-over-ground  problem  with  an  exponentially 
varying  air  density 


at(z)  *  ot(o)e~z/sh  (258) 

atr  "  °tr(°)e"Z/Sh  (259> 

°r  *  ar^°^e"Z/,S^  (260) 

where  at(o),  a®r(o)  are  cross-sections  of  air  at  sea-level. 


z  =  the  height  above  sea-level 
sh  =  atmospheric  scale  height  ~  7km 


In  Eqs  (227)  to  (254)  the  integrals  are  separated  in 
the  space  and  angle  variables.  These  are  double  integrals 
in  space  and  angle.  However,  they  can  be  separated  into 
single  integrals  of  the  y,  x»  P  and  z  variables. 
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APPENDIX  G 


INTEGRALS  FOR  THE 
AIR-OVER- GROUND  PROBLEM 


An  expansion  of  the  air -over -ground  problem  equation, 
Eq  (126),  has  produced  twenty-eight  integral  terms 
(Appendix  F) .  By  a  further  expansion  and  separation  of  the 
integration  variables  thirty-seven  distinct  single  integrals 
are  formed.  These  integrals,  to  include  the  source  inte¬ 
grals  ,  were  numerically  integrated  for  each  combination 
of  the  expansion  subscripts.  They  were  then  stored  in 
a  matrix  and  selected  products  were  used  to  assemble 
the  problem  matrices  and  source  vectors. 

The  thirty-seven  angle  and  spatial  integrals  are 

The  Angle  Integrals 

(261) 

(262) 

(263) 


q cos ( x ) cos  (m  x)cos(nx)dx 

/•  2TT 

Jq  cos(x)|^-|sin(x)cos(mX)|  cos(nX)dX 


•'n 


cos(x)cos(mx)cos (nx)dx 


(264) 


27T 

|jL{sin(x)cos(mX)}  •  |^{sin(x)cos(nx)}dx 

jf  cos  (mx)  •  -|^|sin(x)  -cos(nx)}dx 

(264) 

(265) 

/  cos (mx)cos (nx)dx 

''o 

(266) 

1  cos(x)cos(nx)sin(mx)dx 

JQ 

(267) 

J  cos(nx) *sin(mx)dx 

0 

(268) 

|cos(nx)sin(x)Jsin(mx)dx 

(269) 

f^z  cos  (x )  cos  (mx )  cos  (nx )  dx 

(270) 

/  2cos(x)cos(mx)cos(nx)dx 

(271) 

*v 

■'Z 

r21T 

/  cos(x)cos(mx)cos(nx)dx 

J*y2 

(272) 

fqv?l  m(l,)Pkn(l',du 

(273) 

X1>‘pim<l,>Pkn(l,)dl‘ 

(274) 

r1<1-w2)pkn(M)pJm(P>dv 

(275) 

(276) 


f+1 

-  _ 

C' 

r 

v2  *  lJ*^JLm^^kn^y^dy 

ri 

r 

S_lV*riml  'J>Pta(V»dV 


(277) 

(278) 

(279) 

(280) 


PAtQ(u)  and  Pkn  (y)  are  the  associated  Legendre  polynomials. 


(p)dp 


(281) 

(282) 

(283) 

(284) 

(285) 
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B^r(p) *Bir(P>  pdp 


(286) 


JB 


f0  ez/shBjz(z)Biz(2)dz 


l  ez/sh  |j  {Btz(z)>.Bj2(z)dz 

/i. &.<■>}• 


f 


Biz(z)-Bjz(z)e"Z/Shdz 


H 


f  B.  (z)-B. (z)dz 
Jn  12  J 


The  Source  Integrals 


/fe{vp)}vp)dp 

rR 

I  B.  (p)-H.(p)dp 
^0  Jr  J 

Bjr  (P)Hj (p)pdp 


(287) 

(288) 

(289) 

(290) 

(291) 

(292) 

(293) 

(294) 
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(295) 


-H 


f"  B .  (z)H. (z)ez/shdz 
•'n  Jz  1 


/h{» jz<4iU)e!,8hdI 


^  Bj2(2).Hl(z)dZ 


where 


B(z)  =  cubic  polynomial  z-spline 
B(p)  =  cubic  polynomial  p- spline 
sh  =  atmospheric  scale  height 
R  =  outer  radius  of  the  problem  cylinder 
H  ■  problem  cylinder  height 

H(z)  and  H(p)  are  the  source  interpolating  functions 
Lagrange  polynomials). 


(296) 


(297) 


(linear 
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APPENDIX  H 


AN  EXPANSION  OF  THE  FIRST  SCATTER 
SOURCE  IN  LEGENDRE  POLYNOMIALS 


as 


In  Chapter  VIII  the  first  scatter  source  was  defined 


S(p,z,Q)  =  (298) 


where 

^(r  ,z,6’ )  -  direct  fluence  of  Chapter  VIII,  Eq  (146) 
os  (.z  )  =  scattering  cross-section 

The  usual  legendre  polynomial  cross-section  expansion 
will  now  be  carried  out.  Also  the  even  and  odd  parity 
first  scatter  source  expressions  of  Chapter  VIII  will  be 

c 

derived.  Expanding  a  in  legendre  polynomials  and  using 
the  addition  theorem 


os(z,$-6') 


(z)P£m(y  '  )JWW  )cos  <m<X-x’> ) 


£=0  m*=0 


(299) 


where  m*  means  that  all  terms  with  an  m  ■  0  subscript  must 
be  divided  by  two,  and 


oJ( z)  *  a®(0)e‘2/sh 


(300) 


From  Figure  8  and  Figure  4  it  is  apparent  that 
y'  =  yd  and  x'  *  0,  therefore 

S(p ,z ,fi)  =  2  <fc,(p,z,fi')e~z/8h  * 


L  i 


(301) 


jc=0m -0 


and  using  the  relationships  of  Eqs  (108)  and  (110b)  and  a 
little  algebra,  the  even  and  odd  parity  first  scatter 
sources  can  be  written  as 

Su(p,z,fi)  =  %[S(p,z,H)  -  S(p,zt-fi)]  =  4>.(p,z,fi')e"z/sh  * 


L  i 

EE 

£«0tn*=( 


(302) 
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